Index: /trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py
===================================================================
--- /trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py	(revision 9251)
+++ /trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py	(revision 9252)
@@ -916,4 +916,50 @@
         return allBreakLines
     
+    #########################################
+    def readListOfRiverWalls(rwfileList):
+        """
+            Take a list with the names of riverwall input files [should be comma-separated x,y,elevation files]
+
+            The input file can optionally have a first line defining some hydraulic parameters. A valid example 
+            is
+
+            Qfactor: 1.5, s1: 0.94 
+            200., 300., 0.5
+            300., 400., 0.7
+            ....and so on..
+        
+            Read their coordinates into a dict with their names, read for use by ANUGA
+    
+            INPUT: rwfileList -- a list of riverwall filenames [e.g. from glob.glob('GIS/RiverWalls/*.csv')]
+    
+            OUTPUT: 
+                dictionary with riverwalls [filenames are keys] AND
+                dictionary with hydraulic parameters [filenames are keys]
+        """
+        import numpy
+    
+        allRiverWalls={}
+        allRiverWallPar={}
+        for rwfile in rwfileList:
+            f=open(rwfile)
+            firstLine=f.readline()
+            # If the top line has any letters, assume it is a hydraulic variables line
+            hasLetters=any(c.isalpha() for c in firstLine)
+            if(not hasLetters):
+                allRiverWalls[rwfile]=numpy.genfromtxt(rwfile,delimiter=",").tolist()
+                allRiverWallPar[rwfile]={}
+            else:
+                # Get the wall geometry
+                allRiverWalls[rwfile]=numpy.genfromtxt(rwfile,delimiter=",",skip_header=1).tolist()
+                # Get the hydraulic par
+                firstLine=firstLine.replace(' ', '') # No whitespace
+                wallPar=firstLine.split(',')
+                allRiverWallPar[rwfile]={}
+                for wp in wallPar:
+                    keyNameValue=wp.split(':')
+                    allRiverWallPar[rwfile][keyNameValue[0]]=float(keyNameValue[1])
+        
+        return allRiverWalls, allRiverWallPar
+    
     ############################################################################
     def polygon_from_matching_breaklines(pattern,breakLinesIn, reverse2nd=None):
Index: /trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py
===================================================================
--- /trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py	(revision 9251)
+++ /trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py	(revision 9252)
@@ -10,15 +10,5 @@
 import numpy
 from anuga.structures.inlet_operator import Inlet_operator
-
 from anuga import create_domain_from_regions
-#from anuga import *
-#from swb_domain import domain
-#from anuga import *
-#from balanced_basic import *
-#from balanced_dev import *
-#from bal_and import *
-#from anuga_tsunami import *
-#from balanced_basic.swb2_domain import *
-#from balanced_basic.swb2_domain import Domain
 #------------------------------------------------------------------------------
 # Useful parameters for controlling this case
@@ -30,8 +20,8 @@
 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
+chan_width = 10. # 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)
+l0 = 2.000 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)
 
 assert chan_width < floodplain_width, \
@@ -56,39 +46,48 @@
                      ]
 
-# 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 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. 
 
+breakLines={'n1':[[floodplain_width/2.- chan_width/2., floodplain_length],
+                 [floodplain_width/2. - chan_width/2.,0.]],
+            'n2':[[floodplain_width/2.+ chan_width/2., floodplain_length],
+                  [floodplain_width/2. + chan_width/2.,0.]],
+            # These ones are inside the channel 
+            'n3':[[floodplain_width/2.- chan_width/2.+bankwidth, floodplain_length],
+                 [floodplain_width/2. - chan_width/2.+bankwidth,0.]],
+            'n4':[[floodplain_width/2.+ chan_width/2.-bankwidth, floodplain_length],
+                  [floodplain_width/2. + chan_width/2-bankwidth,0.]]
+            }
+
+regionPtAreas=[[0.01, 0.01, 0.5*l0*l0],
+               [floodplain_width/2.-chan_width/2.+0.01, 0.01, 0.5*l0*l0*0.25],
+               [floodplain_width/2.-chan_width/2.+bankwidth+0.01, 0.01, 0.5*l0*l0],
+               [floodplain_width/2.+chan_width/2.-bankwidth+0.01, 0.01, 0.5*l0*l0*0.25],
+               [floodplain_width/2.+chan_width/2.+0.01, 0.01, 0.5*l0*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)
+anuga.create_mesh_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,
+                         filename = 'channel_floodplain1.msh',
+                         breaklines=breakLines.values(),
+                         regionPtArea=regionPtAreas,
+                         use_cache=False,
+                         verbose=True)
 
+domain=anuga.create_domain_from_file('channel_floodplain1.msh')
 
 domain.set_name('channel_floodplain1') # Output name
-#domain.set_store_vertices_uniquely(True)
+domain.set_store_vertices_uniquely(True)
 domain.set_flow_algorithm('DE1')
-#domain.use_edge_limiter=False
-#domain.extrapolate_velocity_second_order=False
 #------------------------------------------------------------------------------
 #
@@ -100,8 +99,4 @@
 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.))*\
@@ -131,18 +126,7 @@
     return -y*floodplain_slope -chan_bankfull_depth + chan_initial_depth 
 
-domain.set_quantity('elevation', topography) # Use function for elevation
+domain.set_quantity('elevation', topography,location='centroids') # Use function for elevation
 domain.set_quantity('friction', man_n) # Constant friction
 domain.set_quantity('stage', stagetopo) # Use function for stage
-
-#def steady_uh_calc():
-#    # Calculate a boundary condition potentially suitable for steady state flows.
-#    uh = (chan_initial_depth**(4./3.)*man_n**(-2.)*floodplain_slope)**0.5\
-#            *chan_initial_depth  # (u) * h
-#    return uh
-#
-#momin = steady_uh_calc()
-#velin = momin/chan_initial_depth
-#
-#print 'momentum inflow is = ' + str(momin)
 
 # Define inlet operator 
@@ -152,11 +136,9 @@
               [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
+    Qin=4.5
     
     Inlet_operator(domain, line1, Qin)
     
-    print 'Discharge in = ', Qin #,'Velocity at inlet should be ', Qin/(chan_width*chan_initial_depth), \
-            #'for rectangular channels, and approximately that for trapezoidal ones'
+    print 'Discharge in = ', Qin 
 #------------------------------------------------------------------------------
 #
@@ -165,41 +147,6 @@
 #------------------------------------------------------------------------------
 
-# We wish to define one inflow boundary at the upstream head of the channel,
-# one at the outflow of the channel (e.g. transmissive)
-# one at the floodplain edges, excluding the downstream floodplain edge
-#(e.g. reflective)  and one at the downstream floodplain edge (e.g. transmissive)
-
-#Bin_sub = anuga.shallow_water.boundaries.Dirichlet_discharge_boundary(domain,\
-#    chan_initial_depth - chan_bankfull_depth , 0.) # An inflow for subcritical
-#flow
-#Bin_sup = anuga.shallow_water.boundaries.Dirichlet_discharge_boundary(domain, \
-#        chan_initial_depth - chan_bankfull_depth , steady_uh_calc()) # An inflow  
-#                                                                     # for 
-#                                                                     #supercritical 
-#                                                                     # flow
-#Bin_sub = anuga.Dirichlet_boundary([chan_initial_depth - chan_bankfull_depth , \
-#    0., 0.]) # An inflow for subcritical flow -- maybe should use transmissive
-             # boundary instead
-#Bin_sup = anuga.Dirichlet_boundary([chan_initial_depth - chan_bankfull_depth ,\
-#        0., steady_uh_calc()]) # An inflow for supercritical flow
-
-# Trial a novel boundary, which variously imposes 2 or 3 conserved variables
-# depending on the local characteristic speed
-#Bin_trial = anuga.shallow_water.boundaries.\
-#        Dirichlet_inflow_boundary_sub_or_super_critical(domain, \
-#                stage0=chan_initial_depth - chan_bankfull_depth, wh0 = momin,\
-#                href = chan_initial_depth)
-
-# Trial a set_stage_transmissive_momentum inflow boundary
-#def inflow_stage_boundary(t):
-#    return chan_initial_depth - chan_bankfull_depth
-
-# transmissive_momentum_set_stage
-#Bin_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain=domain, function = inflow_stage_boundary) 
-
 Br = anuga.Reflective_boundary(domain) # Solid reflective wall
 Bt = anuga.Transmissive_boundary(domain) # Transmissive boundary
-
-
 Bout_sub = anuga.Dirichlet_boundary( \
         [-floodplain_length*floodplain_slope - chan_bankfull_depth + \
@@ -210,7 +157,4 @@
             + 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 really prevent us reaching steady, uniform flow.
 Bout_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain, function = outflow_stage_boundary) 
 
@@ -224,10 +168,4 @@
                      'chan_in': Br})
 
-
-# To monitor the boundary conditions, let's find centroid points representing sub and super critical flows.
-#inpt = numpy.argmin( ( abs(domain.centroid_coordinates[:,0] - floodplain_width/2.) + \
-#        abs(domain.centroid_coordinates[:,1] - 0.0)  ) )
-#outpt = numpy.argmin( ( abs(domain.centroid_coordinates[:, 0] - floodplain_width/2.) +\
-#        abs(domain.centroid_coordinates[:, 1] - floodplain_length)  ) )
 
 #------------------------------------------------------------------------------
@@ -257,28 +195,3 @@
         print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
 
-    #vv = domain.get_flow_through_cross_section
-    #print domain.quantities['ymomentum'].get_integral()
-    #print 'Qin = ', Qin
-    #print domain.quantities['stage'].centroid_values[inpt]
-    #print domain.quantities['xmomentum'].centroid_values[inpt]
-    #print domain.quantities['ymomentum'].centroid_values[inpt]
 
-    # Alter inflow boundary for sub and super critical flow
-    # FIXME: Varying the inflow boundary conditions does not work well with the modified ANUGA 
-    # (but works great with the standard ANUGA). Especially the supercritical inflow seems bad.
-    # With the modified ANUGA, they oscillate (I wonder if the results are better with a finer grid?) 
-    #h= domain.quantities['stage'].centroid_values[inpt] - domain.quantities['elevation'].centroid_values[inpt]
-    #uh=domain.quantities['xmomentum'].centroid_values[inpt]  
-    #if( domain.g*h < (uh/h)**2):
-    #    domain.set_boundary({'chan_in': Bin_sup})
-    #else:
-    #    domain.set_boundary({'chan_in': Bin_sub})
-
-    # Alter outflow boundary for sub and super critical flow
-    #h=domain.quantities['stage'].centroid_values[outpt] - domain.quantities['elevation'].centroid_values[outpt]
-    #uh=domain.quantities['xmomentum'].centroid_values[outpt]  
-    #if( domain.g*h < (uh/h)**2):
-    #    domain.set_boundary({'chan_out': Bt})
-    #else:
-    #    domain.set_boundary({'chan_out': Bout_sub})
-
Index: /trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py
===================================================================
--- /trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py	(revision 9251)
+++ /trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py	(revision 9252)
@@ -1,33 +1,13 @@
-#import util
-#from bal_and import plot_utils as util
 from anuga.utilities import plot_utils as util
 from matplotlib import pyplot as pyplot
-#import pylab
 
 # Time-index to plot outputs from
-#index=40
-index=570
-#index=225
+index=900
 
-#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_bal_dev_lowvisc.sww', 0.01)
 p2 = util.get_output('channel_floodplain1.sww', 0.001)
 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])
+v = (p.x>5.0)*(p.x<9.0)
 
 # Numerical results along a central channel 'slice'
@@ -40,8 +20,5 @@
 ##########################################################################
 
-#Qin=6.6339  # Inflow discharge
-#Qin=5.2
-Qin=4.6932
-#Qin=4.693286 # Inflow discharge
+Qin=4.5 # Inflow discharge
 slp=1./300. # Floodplain slope (= water slope for steady uniform flow)
 man_n=0.03  # Manning's n
