Index: /trunk/anuga_core/documentation/user_manual/anuga_user_manual.tex
===================================================================
--- /trunk/anuga_core/documentation/user_manual/anuga_user_manual.tex	(revision 8426)
+++ /trunk/anuga_core/documentation/user_manual/anuga_user_manual.tex	(revision 8427)
@@ -1147,6 +1147,7 @@
 if project.scenario == 'fixed_wave':
     # Huge 50m wave starting after 60 seconds and lasting 1 hour.
-    Bw = anuga.Time_boundary(domain=domain,
-                       function=lambda t: [(60<t<3660)*50, 0, 0])
+    Bw = anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(
+                        domain=domain, 
+                        function=lambda t: [(60<t<3660)*50, 0, 0])
     domain.set_boundary({'ocean_east': Bw,
                          'bottom': Bs,
Index: /trunk/anuga_core/documentation/user_manual/demos/cairns/runcairns.py
===================================================================
--- /trunk/anuga_core/documentation/user_manual/demos/cairns/runcairns.py	(revision 8426)
+++ /trunk/anuga_core/documentation/user_manual/demos/cairns/runcairns.py	(revision 8427)
@@ -101,6 +101,8 @@
 if project.scenario == 'fixed_wave':
     # Huge 50m wave starting after 60 seconds and lasting 1 hour.
-    Bw = anuga.Time_boundary(domain=domain,
-                       function=lambda t: [(60<t<3660)*50, 0, 0])
+    Bw = anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(
+                        domain=domain, 
+                        function=lambda t: [(60<t<3660)*50, 0, 0])
+
     domain.set_boundary({'ocean_east': Bw,
                          'bottom': Bs,
Index: /trunk/anuga_core/documentation/user_manual/version.tex
===================================================================
--- /trunk/anuga_core/documentation/user_manual/version.tex	(revision 8426)
+++ /trunk/anuga_core/documentation/user_manual/version.tex	(revision 8427)
@@ -7,3 +7,3 @@
 % release version; this is used to define the
 % \version macro
-\release{1.2.0}
+\release{1.2.1}
Index: /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/channel_floodplain1.py
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/channel_floodplain1.py	(revision 8427)
+++ /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/channel_floodplain1.py	(revision 8427)
@@ -0,0 +1,204 @@
+"""Simple water flow example using ANUGA
+Water flowing down a channel with a floodplain
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+# Import standard shallow water domain and standard boundaries.
+import anuga
+import numpy
+from anuga.structures.inlet_operator import Inlet_operator
+from anuga import *
+#from swb_domain import domain
+#from anuga import *
+#from balanced_basic import *
+#from balanced_dev import *
+#from balanced_basic.swb2_domain import *
+#from balanced_basic.swb2_domain import Domain
+#------------------------------------------------------------------------------
+# Useful parameters for controlling this case
+#------------------------------------------------------------------------------
+
+floodplain_length = 1000.0 # Model domain length
+floodplain_width = 14.0 # Model domain width
+floodplain_slope = 1./300.
+chan_initial_depth = 0.65 # Initial depth of water in the channel
+chan_bankfull_depth = 1.0 # Bankfull depth of the channel
+chan_width = 10.0 # Bankfull width of the channel
+bankwidth = 2. # Width of the bank regions -- note that these protrude into the channel
+man_n=0.03 # Manning's n
+l0 = 1.000 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)
+
+assert chan_width < floodplain_width, \
+        ' ERROR: Channel width is greater than floodplain width'
+
+assert bankwidth < chan_width/2., \
+        'ERROR: The bank width must be less than half the channel width'
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+
+# Define boundary polygon -- in this case clockwise around the proposed boundary
+boundary_polygon = [ [0.,0.], 
+                     [0., floodplain_length], 
+                     [floodplain_width/2. - chan_width/2., floodplain_length], 
+                     [floodplain_width/2. + chan_width/2., floodplain_length], 
+                     [floodplain_width, floodplain_length], 
+                     [floodplain_width, 0.], 
+                     [floodplain_width/2. + chan_width/2., 0.], 
+                     [floodplain_width/2. - chan_width/2., 0.] 
+                     ]
+
+# Define channel polygon, where we can optionally refine the resolution. 
+# Note that we set this a distance l0 inside the boundary polygon, so the polygons
+# do not overlap. 
+channel_polygon = [ [floodplain_width/2. - chan_width/2., +l0],
+                    [floodplain_width/2. - chan_width/2., floodplain_length-l0],
+                    [floodplain_width/2. + chan_width/2., floodplain_length-l0],
+                    [floodplain_width/2. + chan_width/2., +l0]
+                    ]
+
+
+# Define domain with appropriate boundary conditions
+domain = create_domain_from_regions( boundary_polygon, 
+                                   boundary_tags={'left': [0],
+                                                  'top1': [1],
+                                                  'chan_out': [2],
+                                                  'top2': [3],
+                                                  'right': [4],
+                                                  'bottom1': [5],
+                                                  'chan_in': [6],
+                                                  'bottom2': [7] },
+                                   maximum_triangle_area = 0.5*l0*l0,
+                                   minimum_triangle_angle = 28.0,
+                                   mesh_filename = 'channel_floodplain1.msh',
+                                   interior_regions = [ ],
+                                   #interior_regions = [\
+                                   #    [channel_polygon, 0.5*l0*l0] ],
+                                   use_cache=False,
+                                   verbose=True)
+
+
+domain.set_name('channel_floodplain1') # Output name
+#domain.set_store_vertices_uniquely(True)
+#domain.use_edge_limiter=False
+#domain.extrapolate_velocity_second_order=False
+#------------------------------------------------------------------------------
+#
+# Setup initial conditions
+#
+#------------------------------------------------------------------------------
+
+# Function for topography
+def topography(x, y):
+    # Longitudinally sloping floodplain with channel in centre
+	#return -y*floodplain_slope -chan_bankfull_depth*\
+    #        (x>(floodplain_width/2. - chan_width/2.) )*\
+    #        (x<(floodplain_width/2. + chan_width/2.) ) 
+
+    elev1= -y*floodplain_slope - chan_bankfull_depth*\
+            (x>(floodplain_width/2. - chan_width/2.))*\
+            (x<(floodplain_width/2. + chan_width/2.)) 
+    # Add banks
+    if(bankwidth>0.0):
+        leftbnk = floodplain_width/2. - chan_width/2.
+        rightbnk = floodplain_width/2. + chan_width/2.
+        # Left bank
+        elev2 = elev1 + (chan_bankfull_depth \
+                - chan_bankfull_depth/bankwidth*(x - leftbnk))*\
+                (x>leftbnk)*(x < leftbnk + bankwidth)
+        # Right bank
+        elev2 = elev2 + (chan_bankfull_depth \
+                + chan_bankfull_depth/bankwidth*(x - rightbnk))*\
+                (x>rightbnk-bankwidth)*(x < rightbnk)
+    
+    if(bankwidth==0.0):
+        elev2 = elev1
+    #
+    return elev2
+
+#
+
+#Function for stage
+def stagetopo(x,y):
+    return -y*floodplain_slope -chan_bankfull_depth + chan_initial_depth 
+
+domain.set_quantity('elevation', topography) # Use function for elevation
+domain.set_quantity('friction', man_n) # Constant friction
+domain.set_quantity('stage', stagetopo) # Use function for stage
+#domain.set_minimum_allowed_height(0.01) # 
+
+# Define inlet operator 
+flow_in_yval=0.0
+if True:
+    line1 = [ [floodplain_width/2. - chan_width/2., flow_in_yval],\
+              [floodplain_width/2. + chan_width/2., flow_in_yval] \
+              ]
+    Qin = 0.5*(floodplain_slope*(chan_width*chan_initial_depth)**2.*man_n**(-2.)\
+            *chan_initial_depth**(4./3.) )**0.5
+    
+    Inlet_operator(domain, line1, Qin)
+    
+    print 'Discharge in = ', Qin 
+
+#------------------------------------------------------------------------------
+#
+# Setup boundary conditions
+#
+#------------------------------------------------------------------------------
+
+Br = anuga.Reflective_boundary(domain) # Solid reflective wall
+#Bt = anuga.Transmissive_boundary(domain) # Transmissive boundary
+
+def outflow_stage_boundary(t):
+    return -floodplain_length*floodplain_slope \
+            + chan_initial_depth - chan_bankfull_depth
+
+# Note that the outflow boundary may be slightly incorrect for the trapezoidal channel case, 
+# or incorrect more generally if there are numerical problems. But, in the central regions of
+# the channel, this shouldn't prevent us reaching steady, uniform flow.
+Bout_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain, function = outflow_stage_boundary) 
+
+domain.set_boundary({'left': Br, 
+                     'right': Br, 
+                     'top1': Br, 
+                     'top2': Br, 
+                     'bottom1': Br, 
+                     'bottom2': Br, 
+                     'chan_out': Bout_tmss, 
+                     'chan_in': Br})
+
+# Set up file to record computations of discharge at several points.
+discharge_outfile=open('discharge_outputs.txt', 'w')
+discharge_outfile.write('Time (s)'+","+ 'Discharge@10' + ","+ 'Discharge@700'+","+ 'Discharge@1000' + "\n")
+
+#------------------------------------------------------------------------------
+#
+# Evolve system through time
+#
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep=2.0, finaltime=1200.0):
+    print domain.timestepping_statistics()
+    xx=domain.quantities['ymomentum'].centroid_values
+    dd=(domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values)
+    dd=dd*(dd>0.)
+
+    tmp = xx/(dd+1.0e-06)*(dd>0.0)
+    print tmp.max(), tmp.argmax(), tmp.min(),  tmp.argmin()
+
+    # Compute flow through cross-section -- check that the inflow boundary condition is doing its job
+    # This also provides another useful steady-state check
+    if( numpy.floor(t/100.) == t/100. ):
+        print '#### COMPUTING FLOW THROUGH CROSS-SECTIONS########'
+        s0 = domain.get_flow_through_cross_section([[0., 10.0], [floodplain_width, 10.0]])
+        s1 = domain.get_flow_through_cross_section([[0., floodplain_length-300.0], [floodplain_width, floodplain_length-300.0]])
+        s2 = domain.get_flow_through_cross_section([[0., floodplain_length-1.0], [floodplain_width, floodplain_length-1.0]])
+        
+        print 'Cross sectional flows: ',s0, s1, s2 
+        print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
+        discharge_outfile.write(str(t) + "," +str(s0) + ","+ str(s1) +"," + str(s2) + "\n")        
+
+
+discharge_outfile.close()
Index: /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/plotme.py
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/plotme.py	(revision 8427)
+++ /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/plotme.py	(revision 8427)
@@ -0,0 +1,140 @@
+from anuga.utilities import plot_utils as util
+from matplotlib import pyplot as pyplot
+#import pylab
+
+# Time-index to plot outputs from
+index=600
+
+#p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)
+#p=p2
+#p=util.get_centroids(p2, velocity_extrapolation=True)
+
+p2 = util.get_output('channel_floodplain1.sww')
+p=util.get_centroids(p2, velocity_extrapolation=True)
+
+
+#p2 = util.get_output('channel_floodplain1_test_edge_lim.sww', minimum_allowed_height=0.01)
+#p2 = util.get_output('channel_floodplain1_test.sww', minimum_allowed_height=0.01)
+#p=util.get_centroids(p2, velocity_extrapolation=False)
+
+#p2 = util.get_output('channel_floodplain1_standard.sww', minimum_allowed_height=0.01)
+#p2 = util.get_output('channel_floodplain1_balanced_basic.sww', minimum_allowed_height=0.01)
+#p2 = util.get_output('channel_floodplain1_dev.sww')
+#p=util.get_centroids(p2, velocity_extrapolation=True)
+#p=p2
+v = (p.x>6.0)*(p.x<8.0)
+#util.animate_1D(p.time, p.stage[:, v], p.y[v])
+
+# Numerical results along a central channel 'slice'
+V1 = p.stage[index,v] - p.elev[v]
+V2 = p.yvel[index,v]
+V3 = p.xvel[index,v]
+
+##########################################################################
+# Analytical solution of steady uniform 2D flow in a trapezoidal channel.
+##########################################################################
+
+#Qin=6.6339  # Inflow discharge
+#Qin=5.2
+Qin=4.6932
+#Qin=4.693286 # Inflow discharge
+slp=1./300. # Floodplain slope (= water slope for steady uniform flow)
+man_n=0.03  # Manning's n
+Bcentral=6.0 #Flat bed width of the trapezoidal channel
+alpha=0.5  # Side slope of the trapezoidal banks
+
+k = (slp*(1./man_n)**2)**0.5 # At any point, the analytical solution says U = k*d^(2/3)
+
+# Function to calculate the discharge, given the channel centre depth dc, assuming
+# steady uniform flow
+def discharge_su(dc):
+    if(alpha>0.):
+        out = 2*k*( 3./(8.*alpha)*(dc)**(8./3.)) +Bcentral*k*(dc)**(5./3.)
+    else:
+        out = Bcentral*k*(dc)**(5./3.)
+    
+    return out 
+
+# Function that will be minimized to find the depth associated with discharge Qin
+def minme(dc):
+    q1 = discharge_su(dc)
+    return (q1-Qin)**2.
+
+# Minimise the function mimne, to find the centre depth.
+import scipy.optimize
+dc_analytical = scipy.optimize.fmin(minme, x0=1.0)[0]
+
+
+
+##################################
+# Plots
+##################################
+
+# Analytical solution has U*abs(U)*n^2 / D^(4./3.) = Sf = bed slope
+# Hence, the following two variables should be equal -- I have checked that x velocities are fairly small
+pyplot.clf()
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.y[v], (V2**2)**0.5,'o', label='computed velocity')
+pyplot.plot(p.y[v], V2*0+k*dc_analytical**(2./3.), 'o', label='Analytical velocity')
+pyplot.plot(p.y[v], V1, 'o',label='computed depth')
+pyplot.plot(p.y[v], V1*0. + dc_analytical, 'o', label='analytical depth')
+#pyplot.plot( ( (1./300.)*V1**(4./3.)*(1./0.03)**2.)**0.5,'o', label='Analytical velocity based on computed depth')
+pyplot.title('Mid channel velocities and depths, vs analytical velocities and depths')
+# But in my tests, they are not equal
+pyplot.legend( ('computed velocity', 'Analytical velocity', 'computed depth', 'analytical depth'), loc=4)
+pyplot.savefig('trapz_velocity_downstream_l0_eq_1_EL.png')
+
+# Plot velocity over the cross-section
+v1 = (p.y<500.0)&(p.y>490.0)
+
+pyplot.clf()
+analytical_stage = min(p.elev[v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
+pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
+pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
+
+pyplot.legend( ('computed velocity (m/s)', 'analytical velocity (m/s)', 'bed elevation (m)', 'computed stage (m)', 'analytical_stage (m)') ,loc=10)
+pyplot.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (470 to 500m)' +'\n')
+pyplot.savefig('trapz_velocity_cross_channel_l0_eq_1_EL.png') 
+
+
+# Plot velocity over the cross-section
+v1 = (p.y<800.0)&(p.y>790.0)
+
+pyplot.clf()
+analytical_stage = min(p.elev[v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
+pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
+pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
+
+pyplot.legend( ('computed velocity (m/s)', 'analytical velocity (m/s)', 'bed elevation (m)', 'computed stage (m)', 'analytical_stage (m)') ,loc=10)
+pyplot.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (470 to 500m)' +'\n')
+pyplot.savefig('trapz_velocity_cross_channel_l0_eq_1b_EL.png') 
+
+# Plot velocity over the cross-section
+v1 = (p.y<900.0)&(p.y>890.0)
+
+pyplot.clf()
+analytical_stage = min(p.elev[v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
+pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
+pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
+
+pyplot.legend( ('computed velocity (m/s)', 'analytical velocity (m/s)', 'bed elevation (m)', 'computed stage (m)', 'analytical_stage (m)') , loc=10)
+pyplot.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (870 to 900m)' +'\n')
+pyplot.savefig('trapz_velocity_cross_channel_l0_eq_1c_EL.png') 
Index: /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/produce_results.py
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/produce_results.py	(revision 8427)
+++ /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/produce_results.py	(revision 8427)
@@ -0,0 +1,2 @@
+import channel_floodplain1 
+import plotme
Index: /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/results.tex
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/results.tex	(revision 8427)
+++ /trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/results.tex	(revision 8427)
@@ -0,0 +1,27 @@
+\section{Trapezoidal channel steady uniform flow case}
+This scenario simulates steady, uniform flow in a trapezoidal channel. It tests the ability of the model to compute the correct flow depth, and the correct lateral distribution of flow velocity.
+
+\subsection{Results}
+Figure~\ref{downstream_depth} shows the along-channel depth and velocity (in the deepest part of the channel). Both should be nearly constant away from the boundaries, and agree with the analytical solution. 
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{trapz_velocity_downstream_l0_eq_1_EL.png}
+\caption{Depth and velocity mid-channel, in the downstream direction}
+\label{downstream_depth}
+\end{center}
+\end{figure}
+
+Figure~\ref{xsect_vels} shows the cross-channel velocity profiles at a number of cross-sections. Ideally it should agree with the analytical solution, however, this may be difficult due to numerical diffusion in the cross-channel direction. Irrespective, the velocity profile should be qualitatively correct -- highest velocities should be in the channel centre, with lower velocities towards the banks.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{trapz_velocity_cross_channel_l0_eq_1_EL.png}
+\includegraphics[width=0.9\textwidth]{trapz_velocity_cross_channel_l0_eq_1b_EL.png}
+\includegraphics[width=0.9\textwidth]{trapz_velocity_cross_channel_l0_eq_1c_EL.png}
+\caption{Velocity distribution over a number of cross-sections}
+\label{xsect_vels}
+\end{center}
+\end{figure}
+
+\endinput
Index: /trunk/anuga_core/validation_tests/report.tex
===================================================================
--- /trunk/anuga_core/validation_tests/report.tex	(revision 8426)
+++ /trunk/anuga_core/validation_tests/report.tex	(revision 8427)
@@ -30,4 +30,7 @@
 
 \inputresults{Tests/Analytical/runup_sinusoid}
+
+\inputresults{Tests/Analytical/trapezoidal_channel}
+
 %======================
 \chapter{Experimental Tests}
