Index: /trunk/anuga_validation/Hinwood_2008/README.txt
===================================================================
--- /trunk/anuga_validation/Hinwood_2008/README.txt	(revision 8066)
+++ /trunk/anuga_validation/Hinwood_2008/README.txt	(revision 8067)
@@ -4,5 +4,5 @@
 wave tank simulation.
 
-8 experiments were carried out in the flume tank.  Then were four
+8 experiments were carried out in the flume tank.  There were four
 scenarios, with a repeat for each scenario.
 
@@ -20,5 +20,5 @@
   and the given boundary condition and store the model outputs in ANUGA 
   sww files which can be viewed using one of the viewers,
-  e.g. Swollen.  The output is stored in a directory called output.
+  e.g. anuga-viewewr.  The output is stored in a directory called output.
      
 calc_rmds.py
Index: unk/anuga_validation/automated_validation_tests/flow_tests/test_inflow_using_flowline.py
===================================================================
--- /trunk/anuga_validation/automated_validation_tests/flow_tests/test_inflow_using_flowline.py	(revision 8066)
+++ 	(revision )
@@ -1,157 +1,0 @@
-
-
-"""test_inflow_using_flowline
-Test the ability of a flowline to match inflow above the flowline by
-creating constant inflow onto a circle at the head of a 20m
-wide by 400m long plane dipping at various slopes with a perpendicular flowline and gauge 
-downstream of the inflow and a 45 degree flowlines at 200m downstream
-"""
-
-verbose = True
-import numpy as num
-
-#------------------------------------------------------------------------------
-# Import necessary modules
-#------------------------------------------------------------------------------
-
-# Get the standard ANUGA API.
-import anuga
-
-
-#------------------------------------------------------------------------------
-# Setup computational domain
-#------------------------------------------------------------------------------
-
-number_of_inflows = 2 # Number of inflows on top of each other
-finaltime = 1000.0
-
-length = 400.
-width  = 20.
-dx = dy = 5.0          # Resolution: of grid on both axes
-
-points, vertices, boundary = anuga.rectangular_cross(int(length/dx), \
-                                int(width/dy), len1=length, len2=width)
-
-for mannings_n in [0.012, 0.035, 0.070, 0.150]:
-    # Loop over a range of roughnesses  
-    
-    for slope in [1.0/300, 1.0/150, 1.0/75]:
-        # Loop over a range of bedslopes  
-        
-
-        domain = anuga.Domain(points, vertices, boundary)   
-        domain.set_name('inflow_flowline_test')
-
-
-        #----------------------------------------------------------------------
-        # Setup initial conditions
-        #----------------------------------------------------------------------
-
-        def topography(x, y):
-            z=-x * slope
-            return z
-
-        domain.set_quantity('elevation', topography)
-        domain.set_quantity('friction', mannings_n)
-        domain.set_quantity('stage',
-                            expression='elevation') # Dry initial condition
-
-
-        #----------------------------------------------------------------------
-        # Seup Inflow and compute flow characteristics
-        #----------------------------------------------------------------------
-
-        # Fixed Flowrate onto Area 
-        fixed_inflow = anuga.Inflow(domain,
-                              center=(10.0, 10.0),
-                              radius=5.00,
-                              rate=10.00)   
-        
-        # Stack this flow
-        for i in range(number_of_inflows):
-            domain.forcing_terms.append(fixed_inflow)
-        
-        ref_flow = fixed_inflow.rate*number_of_inflows
-        
-        # Compute normal depth on plane using Mannings equation
-        # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
-        normal_depth=(ref_flow*mannings_n/(slope**0.5*width))**0.6
-        flow_velocity = ref_flow/(width*normal_depth)
-        froude_no = flow_velocity/(9.8*normal_depth)**0.5
-        if verbose:
-            print
-            print 'Slope:', slope, 'Mannings n:', mannings_n, 'Velocity:',flow_velocity,'Froude:',froude_no
-
-
-        #----------------------------------------------------------------------
-        # Setup boundary conditions
-        #----------------------------------------------------------------------
-
-        Br = anuga.Reflective_boundary(domain) # Solid reflective wall
-
-        # Define downstream boundary based on predicted depth
-        def normal_depth_stage_downstream(t):
-            return (-slope*length) + normal_depth
-        
-        Bt = anuga.Transmissive_momentum_set_stage_boundary(domain=domain,
-                                                      function=normal_depth_stage_downstream)
-        
-        domain.set_boundary({'left': Br, 'right': Bt, 'top': Br, 'bottom': Br})
-
-        #----------------------------------------------------------------------
-        # Evolve system through time
-        #----------------------------------------------------------------------
-
-        x=200.0      # Reference location for flowlines and gauges
-        y=10.00
-        
-        for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
-            if verbose :
-                print domain.timestepping_statistics()
-                print domain.volumetric_balance_statistics()                
-
-        anuga.sww2dem('inflow_flowline_test.sww',
-        	name_out='sww2dem_out.asc',
-			quantity='xmomentum',
-			cellsize=0.075,      
-			easting_min=0,
-			easting_max=400,
-			northing_min=0,
-			northing_max=20,        
-			reduction=max, 
-			verbose=True)
-
-
-        #----------------------------------------------------------------------
-        # Compute flow thru flowlines ds of inflow and check against normal depth
-        #----------------------------------------------------------------------
-
-        # Square on flowline at 200m
-        q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]])
-        msg = 'Predicted square on flow was %f, should have been %f' % (q, ref_flow)
-        if verbose:
-            print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
-        assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
-        
-        # 45 degree flowline at 200m
-        q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]])
-        msg = 'Predicted 45 flow was %f, should have been %f' % (q, ref_flow)
-        if verbose:
-            print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
-        assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
-
-
-       # Stage recorder (gauge) in middle of plane at 200m
-        w = domain.get_quantity('stage').get_values(interpolation_points=[[x, y]])[0]
-        z = domain.get_quantity('elevation').get_values(interpolation_points=[[x, y]])[0]
-        domain_depth = w-z
-        
-        
-        msg = 'Predicted depth of flow at 200m was %f, should have been %f' % (normal_depth, domain_depth)
-        msg += 'Mannings n = %s, slope = %s' % (mannings_n, slope)
-        if verbose:
-            print 'Depth at 200m: ANUGA = %f, Mannings = %f' % (domain_depth, normal_depth)
-        assert num.allclose(domain_depth,normal_depth, rtol=1.0e-2), msg
-
-        
-
Index: /trunk/anuga_validation/automated_validation_tests/flow_tests/validate_flow.py
===================================================================
--- /trunk/anuga_validation/automated_validation_tests/flow_tests/validate_flow.py	(revision 8066)
+++ /trunk/anuga_validation/automated_validation_tests/flow_tests/validate_flow.py	(revision 8067)
@@ -25,7 +25,7 @@
         
         #print
-        s = 'test_inflow_using_flowline.py'
+        s = 'run_inflow_test_using_flowline.py'
         #print s
-        res = os.system('python %s > test_inflow_using_flowline.stdout' %s)
+        res = os.system('python %s > run_inflow_test_using_flowline.stdout' %s)
         assert res == 0
 
Index: /trunk/anuga_validation/automated_validation_tests/validate_all.py
===================================================================
--- /trunk/anuga_validation/automated_validation_tests/validate_all.py	(revision 8066)
+++ /trunk/anuga_validation/automated_validation_tests/validate_all.py	(revision 8067)
@@ -34,5 +34,5 @@
 
     for filename in filenames:
-        if filename.startswith('validate') and filename.endswith('.py'):
+        if filename.startswith('validate_') and filename.endswith('.py'):
             #print 'Found %s in %s' %(filename, dirpath)
             validation_dirs_and_files.append((dirpath, filename))            
Index: /trunk/anuga_validation/convergence_study/tilted.py
===================================================================
--- /trunk/anuga_validation/convergence_study/tilted.py	(revision 8066)
+++ /trunk/anuga_validation/convergence_study/tilted.py	(revision 8067)
@@ -11,15 +11,15 @@
 
 import sys
-from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-from anuga.shallow_water import Domain
-from anuga.shallow_water import Reflective_boundary
-from anuga.shallow_water import Dirichlet_boundary
-from anuga.shallow_water import Time_boundary
-from anuga.shallow_water import Transmissive_boundary
-from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
-from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+from anuga import Transmissive_boundary
+from anuga import Transmissive_momentum_set_stage_boundary
+#from anuga.data_manager import start_screen_catcher, copy_code_files
 
 from math import cos
-from Numeric import zeros, Float
+import numpy as num
 from time import localtime, strftime, gmtime
 from os import sep
@@ -36,5 +36,5 @@
 output_dir = model_name + '_'+time
 output_file = model_name
-copy_code_files(output_dir,__file__)
+#copy_code_files(output_dir,__file__)
 #start_screen_catcher(output_dir+sep)
 
@@ -63,6 +63,6 @@
 print domain.get_timestepping_method()
 
-domain.use_edge_limiter = True
-domain.tight_slope_limiters = True
+#domain.use_edge_limiter = True
+#domain.tight_slope_limiters = True
 
 domain.CFL = 1.0
Index: /trunk/anuga_validation/convergence_study/wave.py
===================================================================
--- /trunk/anuga_validation/convergence_study/wave.py	(revision 8066)
+++ /trunk/anuga_validation/convergence_study/wave.py	(revision 8067)
@@ -11,15 +11,15 @@
 
 import sys
-from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-from anuga.shallow_water import Domain
-from anuga.shallow_water import Reflective_boundary
-from anuga.shallow_water import Dirichlet_boundary
-from anuga.shallow_water import Time_boundary
-from anuga.shallow_water import Transmissive_boundary
-from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
-from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+from anuga import Transmissive_boundary
+from anuga import Transmissive_momentum_set_stage_boundary
+#from anuga.data_manager import start_screen_catcher, copy_code_files
 
 from math import cos
-from Numeric import zeros, Float
+import numpy as num
 from time import localtime, strftime, gmtime
 from os import sep
@@ -36,8 +36,8 @@
 output_file = 'wave'
 
-copy_code_files(output_dir,__file__)
+#copy_code_files(output_dir,__file__)
 #start_screen_catcher(output_dir+sep)
 
-interactive_visualisation = False
+interactive_visualisation = True
 
 #------------------------------------------------------------------------------
@@ -60,21 +60,24 @@
 # Setup Algorithm
 #------------------------------------------------------------------------------
-domain.set_timestepping_method('rk2')
+
+domain.set_timestepping_method(2)
 domain.set_default_order(2)
 
 print domain.get_timestepping_method()
 
-domain.use_edge_limiter = True
-domain.tight_slope_limiters = False
-domain.use_centroid_velocities = False
+#domain.use_edge_limiter = True
+#domain.tight_slope_limiters = False
+#domain.use_centroid_velocities = False
+
+domain.set_beta(1.9)
 
 domain.CFL = 1.0
 
-domain.beta_w      = 1.0
-domain.beta_w_dry  = 0.0
-domain.beta_uh     = 1.0
-domain.beta_uh_dry = 0.0
-domain.beta_vh     = 1.0
-domain.beta_vh_dry = 0.0
+#domain.beta_w      = 1.0
+#domain.beta_w_dry  = 0.0
+#domain.beta_uh     = 1.0
+#domain.beta_uh_dry = 0.0
+#domain.beta_vh     = 1.0
+#domain.beta_vh_dry = 0.0
 
 
Index: /trunk/anuga_validation/validation_tests/dam_break/dam_break.py
===================================================================
--- /trunk/anuga_validation/validation_tests/dam_break/dam_break.py	(revision 8067)
+++ /trunk/anuga_validation/validation_tests/dam_break/dam_break.py	(revision 8067)
@@ -0,0 +1,119 @@
+"""Simple water flow example using ANUGA
+
+Water driven up a linear slope and time varying boundary,
+similar to a beach environment
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+import anuga
+
+from math import cos
+from numpy import zeros, float
+from time import localtime, strftime, gmtime
+
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+output_dir = 'dam_break_'+time
+output_file = 'dam_break'
+
+anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 200.
+dy = dx
+L = 100000.
+W = 10*dx
+
+# structured mesh
+points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+domain = anuga.Domain(points, vertices, boundary) 
+
+domain.set_name(output_file)                
+domain.set_datadir(output_dir)  
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+domain.set_timestepping_method('rk2')
+domain.set_default_order(2)
+
+print domain.get_timestepping_method()
+
+domain.use_edge_limiter = True
+domain.tight_slope_limiters = True
+domain.use_centroid_velocities = False
+
+domain.CFL = 1.0
+
+domain.beta_w      = 0.6
+domain.beta_uh     = 0.6
+domain.beta_vh     = 0.6
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+domain.set_quantity('elevation',0.0)
+domain.set_quantity('friction', 0.0)
+
+h0 = 10.0
+h1 = 0.0
+
+def height(x,y):
+    z = zeros(len(x), float)
+    for i in range(len(x)):
+        if x[i]<=50000.0:
+            z[i] = h0
+        else:
+            z[i] = h1
+    return z
+
+
+domain.set_quantity('stage', height)
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+from math import sin, pi, exp
+Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
+Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary 
+Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
+
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bt, 'right': Bt, 'top': Br, 'bottom': Br})
+
+
+#===============================================================================
+from anuga.visualiser import RealtimeVisualiser
+vis = RealtimeVisualiser(domain)
+vis.render_quantity_height("stage", zScale =h0*500, dynamic=True)
+vis.colour_height_quantity('stage', (0.0, 0.5, 1.0))
+vis.start()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 100.0, finaltime = 60*60.):
+    print domain.timestepping_statistics(track_speeds=True)
+    vis.update()
+    
+vis.evolveFinished()
+
Index: /trunk/anuga_validation/validation_tests/periodic_carrier/carrier_wave_runup.py
===================================================================
--- /trunk/anuga_validation/validation_tests/periodic_carrier/carrier_wave_runup.py	(revision 8067)
+++ /trunk/anuga_validation/validation_tests/periodic_carrier/carrier_wave_runup.py	(revision 8067)
@@ -0,0 +1,218 @@
+"""Example of shallow water wave equation analytical solution of the
+one-dimensional Thacker and Greenspan wave run-up treated as a two-dimensional solution.
+
+   Copyright 2004
+   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
+   Geoscience Australia
+
+Specific methods pertaining to the 2D shallow water equation
+are imported from shallow_water
+for use with the generic finite volume framework
+
+Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
+numerical vector named conserved_quantities.
+"""
+
+#-------------------
+# Module imports
+import sys
+from math import sqrt, cos, sin, pi
+from numpy import asarray
+from time import localtime, strftime, gmtime
+
+
+#-------------------
+# Anuga Imports
+#-------------------
+from anuga.shallow_water_balanced.swb_domain import Domain, Transmissive_boundary, Reflective_boundary,\
+     Dirichlet_boundary, Time_boundary
+
+#from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\
+#     Dirichlet_boundary
+
+from anuga.utilities.polygon import inside_polygon, is_inside_triangle
+from anuga.abstract_2d_finite_volumes.mesh_factory import strang_mesh
+
+
+#-------------------
+#Convenience functions
+#-------------------
+def imag(a):
+    return a.imag
+
+def real(a):
+    return a.real
+
+
+#--------------------
+# Domain
+# Strang_domain will search through the file and test to see if there are
+# two or three entries. Two entries are for points and three for triangles.
+
+#points, elements = strang_mesh('Run-up.pt')
+points, elements = strang_mesh('strang_7389.pt')
+domain = Domain(points, elements)
+
+print 'extent ',domain.get_extent()
+
+#----------------
+# Order of scheme
+# Good compromise between
+# limiting and CFL
+#---------------
+domain.set_default_order(2)
+domain.set_timestepping_method(2)
+domain.set_beta(0.7)
+domain.set_CFL(0.6)
+
+
+#-------------------
+#Set a default tagging
+#-------------------
+epsilon = 1.0e-12
+
+for id, face in domain.boundary:
+    domain.boundary[(id,face)] = 'external'
+    if domain.get_vertex_coordinate(id,(face+1)%3)[0] < -200.0+1.0e-10:
+        domain.boundary[(id,face)] = 'left'
+
+#---------------------
+# Provide file name for storing output
+#---------------------
+domain.store = True
+domain.format = 'sww'
+
+time = strftime('%Y%m%d_%H%M%S',localtime())
+output_file= 'carrier_wave_runup_'+time
+
+domain.set_name(output_file)
+
+print "Number of triangles = ", len(domain)
+
+#-----------------------
+#Define the boundary condition
+#-----------------------
+def stage_setup(x,t):
+    vh = 0
+    alpha = 0.1
+    eta = 0.1
+    a = 1.5*sqrt(1.+0.9*eta)
+    l_0 = 200.
+    ii = complex(0,1)
+    g = 9.81
+    v_0 = sqrt(g*l_0*alpha)
+    v1 = 0.
+
+    sigma_max = 100.
+    sigma_min = -100.
+    for j in range (1,50):
+        sigma0 = (sigma_max+sigma_min)/2.
+        lambda_prime = 2./a*(t/sqrt(l_0/alpha/g)+v1)
+        sigma_prime = sigma0/a
+        const = (1.-ii*lambda_prime)**2+sigma_prime**2
+
+        v1 = 8.*eta/a*imag(1./const**(3./2.) \
+            -3./4.*(1.-ii*lambda_prime)/const**(5./2.))
+
+        x1 = -v1**2/2.-a**2*sigma_prime**2/16.+eta \
+            *real(1.-2.*(5./4.-ii*lambda_prime) \
+            /const**(3./2.)+3./2.*(1.-ii*lambda_prime)**2 \
+            /const**(5./2.))
+
+        neta1 = x1 + a*a*sigma_prime**2/16.
+
+        v_star1 = v1*v_0
+        x_star1 = x1*l_0
+        neta_star1 = neta1*alpha*l_0
+        stage = neta_star1
+        z = stage - x_star1*alpha
+        uh = z*v_star1
+
+        if x_star1-x > 0:
+            sigma_max = sigma0
+        else:
+            sigma_min = sigma0
+
+        if abs(abs(sigma0)-100.) < 10:
+
+#     solution does not converge because bed is dry
+            stage = 0.
+            uh = 0.
+            z = 0.
+
+    return [stage, uh, vh]
+
+def boundary_stage(t):
+    x = -200
+    return stage_setup(x,t)
+
+
+#---------------------
+#Initial condition
+#---------------------
+print 'Initial condition'
+t_star1 = 0.0
+slope = -0.1
+
+#Set bed-elevation and friction(None)
+def x_slope(x,y):
+    n = x.shape[0]
+    z = 0*x
+    for i in range(n):
+        z[i] = -slope*x[i]
+    return z
+
+domain.set_quantity('elevation', x_slope)
+
+#Set the water depth
+print 'Initial water depth'
+def stage(x,y):
+    z = x_slope(x,y)
+    n = x.shape[0]
+    w = 0*x
+    for i in range(n):
+        w[i], uh, vh = stage_setup(x[i],t_star1)
+        h = w[i] - z[i]
+        if h < 0:
+            h = 0
+            w[i] = z[i]
+    return w
+
+domain.set_quantity('stage', stage)
+
+
+#####################
+#Set up boundary conditions
+Br = Reflective_boundary(domain)
+Bw = Time_boundary(domain, boundary_stage)
+domain.set_boundary({'left': Bw, 'external': Br})
+
+import time
+visualize = True
+if visualize:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("elevation", zScale=3.0, offset = 0.01, dynamic=False)
+    vis.render_quantity_height("stage", zScale = 3.0, dynamic=True, opacity = 0.6, wireframe=False)
+    #vis.colour_height_quantity('stage', (lambda q:q['stage'], 1.0, 2.0))
+    vis.colour_height_quantity('stage', (0.4, 0.6, 0.4))
+    vis.start()
+    time.sleep(2.0)
+
+#domain.visualise = True
+
+######################
+#Evolution
+
+t0 = time.time()
+
+
+for t in domain.evolve(yieldstep = 1., finaltime = 100):
+    domain.write_time()
+    print boundary_stage(domain.time)
+    if visualize: vis.update()
+    
+if visualize: vis.evolveFinished()    
+
+print 'That took %.2f seconds' %(time.time()-t0)
+
Index: /trunk/anuga_validation/validation_tests/periodic_forcing/wave.py
===================================================================
--- /trunk/anuga_validation/validation_tests/periodic_forcing/wave.py	(revision 8067)
+++ /trunk/anuga_validation/validation_tests/periodic_forcing/wave.py	(revision 8067)
@@ -0,0 +1,131 @@
+"""Simple water flow example using ANUGA
+
+Will Powers example of a simple sinusoidal wave which showed diffusive effects of
+thefirst order and standard second order method. Problem resolved if "rk2" timestepping
+and higher beta = 2 limiter used. Also new edge limiter with rk2 resolves problem 
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import Time_boundary
+from anuga.shallow_water import Transmissive_boundary
+from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
+
+from math import cos
+from Numeric import zeros, Float
+from time import localtime, strftime, gmtime
+from os import sep
+
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+output_dir = 'wave_'+time
+output_file = 'wave'
+
+copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+sep)
+
+interactive_visualisation = False
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1000.
+dy = dx
+L = 100000.
+W = 10*dx
+
+# structured mesh
+points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+domain = Domain(points, vertices, boundary) 
+
+domain.set_name(output_file)                
+domain.set_datadir(output_dir)  
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+domain.set_timestepping_method('rk2')
+domain.set_default_order(2)
+
+print domain.get_timestepping_method()
+
+domain.use_edge_limiter = True
+domain.tight_slope_limiters = False
+domain.use_centroid_velocities = False
+
+domain.CFL = 1.0
+
+domain.beta_w      = 1.0
+domain.beta_w_dry  = 0.0
+domain.beta_uh     = 1.0
+domain.beta_uh_dry = 0.0
+domain.beta_vh     = 1.0
+domain.beta_vh_dry = 0.0
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+domain.set_quantity('elevation',-100.0)
+domain.set_quantity('friction', 0.00)
+domain.set_quantity('stage', 0.0)            
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+from math import sin, pi, exp
+Br = Reflective_boundary(domain)      # Solid reflective wall
+Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
+Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
+amplitude = 1
+Bw = Time_boundary(domain=domain,     # Time dependent boundary  
+## Sine wave
+                  f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
+## Sawtooth?
+#                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
+## Sharp rise, linear fall
+#                   f=lambda t: [(5.0*(-((t-0.)/300.)*(t<300.)-cos((t-300.)*2.*pi*(1./240.))*(t>=300. and t<420.)+(1.-(t-420.)/300.)*(t>=420. and t <720.))), 0.0, 0.0])
+#                   f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
+#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bw, 'right': Bt, 'top': Br, 'bottom': Br})
+
+
+#===============================================================================
+if interactive_visualisation:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("stage", zScale =10000, dynamic=True)
+    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
+    vis.start()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 50.0, finaltime = 60*60.):
+    domain.write_time()
+    if interactive_visualisation:
+        vis.update()
+
+if interactive_visualisation:
+    vis.evolveFinished()
+
Index: /trunk/anuga_validation/validation_tests/yoon_parabolic/yoon_import_mesh.py
===================================================================
--- /trunk/anuga_validation/validation_tests/yoon_parabolic/yoon_import_mesh.py	(revision 8067)
+++ /trunk/anuga_validation/validation_tests/yoon_parabolic/yoon_import_mesh.py	(revision 8067)
@@ -0,0 +1,141 @@
+"""Example of shallow water wave equation analytical solution
+consists of a parabolic profile in a parabolic basin. Analytical
+solutiuon to this problem was derived by Carrier and Greenspan
+and used by Yoon and Chou.
+
+   Copyright 2005
+   Christopher Zoppou, Stephen Roberts, Ole Nielsen
+   ANU, Geoscience Australia
+   
+"""
+
+#---------------
+# Module imports 
+import sys
+
+from anuga.interface import Domain, Dirichlet_boundary
+from math import sqrt, cos, sin, pi
+from anuga.abstract_2d_finite_volume.mesh_factory import strang_mesh
+from anuga.utilities.polygon import inside_polygon
+
+from numpy import asarray
+from anuga.fit_interpolate.interpolate import Interpolation
+
+#-------------------------------
+# Set up the domain of triangles
+# Strang_domain will search
+# through the file and test to
+# see if there are two or three
+# entries. Two entries are for
+# points and three for triangles.
+points, elements = strang_mesh('yoon_circle.pt')
+domain = Domain(points, elements) 
+
+#----------------
+# Order of scheme
+domain.default_order = 2
+
+domain.smooth = True
+
+#-------------------------------------
+# Provide file name for storing output
+domain.store = False
+domain.format = 'sww'
+domain.set_name('yoon_mesh_second_order.2')
+print 'Number of triangles = ', len(domain)
+
+#----------------------------------------------------------
+# Decide which quantities are to be stored at each timestep
+domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+#------------------------------------------
+# Reduction operation for get_vertex_values
+from anuga.pyvolution.util import mean
+domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
+
+
+#------------------
+# Initial condition
+print 'Initial condition'
+t = 0.0
+D0 = 1.
+L = 2500.
+R0 = 2000.
+g = 9.81
+
+A = (L**4 - R0**4)/(L**4 + R0**4)
+omega = 2./L*sqrt(2.*g*D0)
+T = pi/omega
+
+#------------------
+# Set bed elevation
+def x_slope(x,y):
+    n = x.shape[0]
+    z = 0*x
+    for i in range(n):
+        r = sqrt(x[i]*x[i] + y[i]*y[i])
+        z[i] = -D0*(1.-r*r/L/L)
+    return z
+domain.set_quantity('elevation', x_slope)
+
+#----------------------------
+# Set the initial water level
+def level(x,y):
+    z = x_slope(x,y)
+    n = x.shape[0]
+    h = 0*x
+    for i in range(n):
+        r = sqrt(x[i]*x[i] + y[i]*y[i])
+        h[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
+                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
+    if h[i] < z[i]:
+        h[i] = z[i]
+    return h   
+domain.set_quantity('stage', level)
+
+
+#---------
+# Boundary
+print 'Boundary conditions'
+domain.set_boundary({'exterior': Dirichlet_boundary([0.0, 0.0, 0.0])})
+
+#---------------------------------------------
+# Find triangle that contains the point points
+# and print to file
+points = [0.,0.]
+for n in range(len(domain.triangles)):
+    t = domain.triangles[n]
+    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+
+    if inside_polygon(points,tri):
+        print 'Point is within triangle with vertices '+'%s'%tri
+        n_point = n
+t = domain.triangles[n_point]
+tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+
+filename=domain.get_name()
+file = open(filename,'w')
+    
+#----------
+# Evolution
+import time
+t0 = time.time()
+for t in domain.evolve(yieldstep = 20.0, finaltime = 3000 ):
+    domain.write_time()
+
+    tri_array = asarray(tri)
+    t_array = asarray([[0,1,2]])
+    interp = Interpolation(tri_array,t_array,[points])
+
+    stage = domain.get_quantity('stage').get_values()[n_point]
+    xmomentum = domain.get_quantity('xmomentum').get_values()[n_point]
+    ymomentum = domain.get_quantity('ymomentum').get_values()[n_point]
+
+    interp_stage = interp.interpolate([[stage[0]],[stage[1]],[stage[2]]])
+    interp_xmomentum = interp.interpolate([[xmomentum[0]],[xmomentum[1]],[xmomentum[2]]])
+    interp_ymomentum = interp.interpolate([[ymomentum[0]],[ymomentum[1]],[ymomentum[2]]])
+    
+    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,interp_stage[0][0],interp_xmomentum[0][0],interp_ymomentum[0][0]) )
+
+file.close()    
+print 'That took %.2f seconds' %(time.time()-t0)
Index: /trunk/anuga_validation/validation_tests/yoon_parabolic/yoon_parabolic_basin_cross_mesh.py
===================================================================
--- /trunk/anuga_validation/validation_tests/yoon_parabolic/yoon_parabolic_basin_cross_mesh.py	(revision 8067)
+++ /trunk/anuga_validation/validation_tests/yoon_parabolic/yoon_parabolic_basin_cross_mesh.py	(revision 8067)
@@ -0,0 +1,177 @@
+"""Example of shallow water wave equation analytical solution
+consists of a parabolic profile in a parabolic basin. Analytical
+solutiuon to this problem was derived by Thacker
+and used by Yoon and Chou.
+
+   Copyright 2005
+   Christopher Zoppou, Stephen Roberts, ANU, Geoscience Australia
+
+"""
+
+#---------------
+# Module imports
+from anuga.shallow_water_balanced.swb_domain import Domain, Transmissive_boundary, Reflective_boundary,\
+     Dirichlet_boundary
+
+#from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\
+#     Dirichlet_boundary
+from math import sqrt, cos, sin, pi
+from anuga.interface import rectangular_cross
+from anuga.utilities.polygon import inside_polygon, is_inside_triangle
+from numpy import asarray
+
+
+#-------------------------------
+# Domain
+n = 200
+m = 200
+lenx = 8000.0
+leny = 8000.0
+origin = (-4000.0, -4000.0)
+
+points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
+domain = Domain(points, elements, boundary)
+
+#----------------
+# Order of scheme
+# Good compromise between
+# limiting and CFL
+#---------------
+domain.set_default_order(2)
+domain.set_timestepping_method(2)
+domain.set_beta(0.7)
+domain.set_CFL(0.6)
+
+domain.smooth = True
+
+#-------------------------------------
+# Provide file name for storing output
+domain.store = False
+domain.format = 'sww'
+domain.set_name('yoon_mesh_second_order_cross')
+print 'Number of triangles = ', len(domain)
+
+#----------------------------------------------------------
+# Decide which quantities are to be stored at each timestep
+domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+#------------------------------------------
+# Reduction operation for get_vertex_values
+from anuga.utilities.numerical_tools import mean
+domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
+
+#------------------
+# Initial condition
+print 'Initial condition'
+t = 0.0
+D0 = 1.
+L = 2500.
+R0 = 2000.
+g = 9.81
+
+A = (L**4 - R0**4)/(L**4 + R0**4)
+omega = 2./L*sqrt(2.*g*D0)
+T = pi/omega
+
+#------------------
+# Set bed elevation
+def x_slope(x,y):
+    n = x.shape[0]
+    z = 0*x
+    for i in range(n):
+        r = sqrt(x[i]*x[i] + y[i]*y[i])
+        z[i] = -D0*(1.-r*r/L/L)
+    return z
+domain.set_quantity('elevation', x_slope)
+
+#----------------------------
+# Set the initial water level
+def level(x,y):
+    z = x_slope(x,y)
+    n = x.shape[0]
+    h = 0*x
+    for i in range(n):
+        r = sqrt(x[i]*x[i] + y[i]*y[i])
+        h[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
+                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
+    if h[i] < z[i]:
+        h[i] = z[i]
+    return h
+domain.set_quantity('stage', level)
+
+#---------
+# Boundary
+print 'Boundary conditions'
+R = Reflective_boundary(domain)
+T = Transmissive_boundary(domain)
+D = Dirichlet_boundary([0.0, 0.0, 0.0])
+domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D})
+
+#---------------------------------------------
+# Find triangle that contains the point points
+# and print to file
+
+
+
+points = (0.0, 0.0)
+for n in range(len(domain.triangles)):
+    tri = domain.get_vertex_coordinates(n)
+
+    if is_inside_triangle(points,tri):
+        #print 'Point is within triangle with vertices '+'%s'%tri
+        n_point = n
+
+print 'n_point = ',n_point
+t = domain.triangles[n_point]
+print 't = ', t
+tri = domain.get_vertex_coordinates(n)
+
+filename=domain.get_name()
+file = open(filename,'w')
+
+#----------
+# Evolution
+import time
+t0 = time.time()
+
+time_array = []
+stage_array = []
+Stage     = domain.quantities['stage']
+Xmomentum = domain.quantities['xmomentum']
+Ymomentum = domain.quantities['ymomentum']
+
+for t in domain.evolve(yieldstep = 20.0, finaltime = 17700.0 ):
+    domain.write_time()
+
+    #tri_array = asarray(tri)
+    #t_array = asarray([[0,1,2]])
+    #interp = Interpolation(tri_array,t_array,[points])
+
+
+    stage     = Stage.get_values(location='centroids',indices=[n_point])
+    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])
+    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])
+    #print '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum)
+    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )
+
+    time_array.append(t)
+    stage_array.append(stage)
+
+file.close()
+print 'That took %.2f seconds' %(time.time()-t0)
+
+
+from pylab import *
+ion()
+hold(False)
+plot(time_array, stage_array, 'r.-')
+#title('Gauge %s' %name)
+xlabel('time(s)')
+ylabel('stage (m)')
+#legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
+#savefig(name, dpi = 300)
+
+#raw_input('Next')
+show()
+
+
