Index: /trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py
===================================================================
--- /trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py	(revision 8382)
+++ /trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py	(revision 8383)
@@ -93,5 +93,5 @@
         # Because gravity is treated within the flux function, 
         # we remove it from the forcing terms.
-        self.forcing_terms.remove(gravity)
+        #self.forcing_terms.remove(gravity)
 
         # We need the edge_coordinates for the extrapolation
Index: /trunk/anuga_work/development/gareth/tests/urban_flow/ideal_urban.py
===================================================================
--- /trunk/anuga_work/development/gareth/tests/urban_flow/ideal_urban.py	(revision 8383)
+++ /trunk/anuga_work/development/gareth/tests/urban_flow/ideal_urban.py	(revision 8383)
@@ -0,0 +1,168 @@
+"""Idealised urban flow example using ANUGA
+Water flowing over a floodplain with buildings 
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+# Import standard shallow water domain and standard boundaries.
+import anuga
+import numpy
+#from anuga.structures.inlet_operator import Inlet_operator
+
+from balanced_dev import *
+from balanced_dev import create_domain_from_regions as create_domain_from_regions
+#------------------------------------------------------------------------------
+# Useful parameters for controlling this case
+#------------------------------------------------------------------------------
+
+floodplain_length = 35.80 # Model domain length
+floodplain_width = 3.60 # Model domain width
+floodplain_slope = 0.0 #1./3000.
+#chan_initial_depth = 0.8 # 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.01 # Manning's n
+l0 = 0.05 #0.2501 # 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_width], 
+                     [floodplain_length, floodplain_width], 
+                     [floodplain_length, 0.], 
+                     ]
+
+
+# Define domain with appropriate boundary conditions
+domain = create_domain_from_regions( boundary_polygon, 
+                                   boundary_tags={'left': [0],
+                                                  'top': [1],
+                                                  'right': [2],
+                                                  'bottom': [3]},
+                                   maximum_triangle_area = 1.0*l0*l0,
+                                   minimum_triangle_angle = 28,
+                                   mesh_filename = 'urban_flow1.msh',
+                                   use_cache=False,
+                                   verbose=True)
+
+
+domain.set_name('urban_flow0p05') # Output name
+
+#------------------------------------------------------------------------------
+#
+# Setup initial conditions
+#
+#------------------------------------------------------------------------------
+
+## Snippet of code to test the topography function
+## We define x,y, then look at z(x,y) to see if it is correct (I use
+## pyplot.scatter for this)
+## Predefine x,y to be long enough
+#increm=0.05
+#x = numpy.ones(int(((floodplain_length)/increm+1.0)*((floodplain_width)/increm+1.0) ))
+#y = numpy.ones(int(((floodplain_length)/increm+1.0)*((floodplain_width)/increm+1.0) ))
+## Fill out x,y
+#counter=0
+#for i in numpy.arange(0.0,floodplain_length, increm):
+#   for j in numpy.arange(0.0, floodplain_width, increm): 
+#       x[counter]=i
+#       y[counter]=j
+#       counter=counter+1
+#
+#z = topography(x,y)
+
+# Function for topography, based on S. Soares Frazao, Noel and Zech,
+# Experiments of dam break flow in the presence of obstables.
+# J. Hydraulic Research
+def topography(x, y):
+    #import pdb
+    #pdb.set_trace()
+    # Define 'bare-earth' topography
+    structure_height=0.8
+    elev1= 0.0*y #-y*floodplain_slope 
+    # Define constriction
+    elev2 = (x>6.9)*(x<7.7)*(y<1.3)*structure_height + \
+            (x>6.9)*(x<7.7)*(3.6-y<1.3)*structure_height
+    # Define building.
+    # Do this by using rotated coordinates.
+    alpha = 64.0/180.0*numpy.pi
+    short_side = 0.4
+    long_side=0.8
+    # The lower left building corner is at (11.10, 1.75)
+    ll = (11.10, 1.75)
+    # The lower right building corner is at (11.10, 1.75) + 0.4*(numpy.sin(alpha), -numpy.cos(alpha))
+    lr = (ll[0] + short_side*numpy.sin(alpha), ll[1] -short_side*numpy.cos(alpha))
+    # The upper left building corner is at (11.10, 1.75) + 0.8*(numpy.cos(64./180.), numpy.sin(64./180.)) 
+    ul = (ll[0] + long_side*numpy.cos(alpha), ll[1] + long_side*numpy.sin(alpha))
+    c1 = y-numpy.tan(alpha)*x # Coordinate parallel to the long side of the building
+    c2 = y +1.0/(numpy.tan(alpha))*x # Coordinate parallel to the short side of the building
+    elev3 = (c1 < (ll[1] - numpy.tan(alpha)*ll[0]) )*\
+            (c1 > (lr[1] - numpy.tan(alpha)*lr[0]) )*\
+            (c2 < (ul[1] + 1.0/numpy.tan(alpha)*ul[0]) )*\
+            (c2 > (ll[1] + 1.0/numpy.tan(alpha)*ll[0]) )*\
+            structure_height
+    return elev1 + elev2 + elev3
+
+#Function for stage
+def stagetopo(x,y):
+    return  (x<7.7)*0.39 + 0.01
+
+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
+
+
+#------------------------------------------------------------------------------
+#
+# Setup boundary conditions
+#
+#------------------------------------------------------------------------------
+
+# 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) 
+
+# Standard boundary conditions
+Br = anuga.Reflective_boundary(domain) # Solid reflective wall
+Bt = anuga.Transmissive_boundary(domain) # Transmissive boundary
+# Outflow boundary conditions
+def outflow_stage_boundary(t):
+    return -floodplain_length*floodplain_slope \
+            + chan_initial_depth - chan_bankfull_depth
+Bout_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain, function = outflow_stage_boundary) 
+
+# Define boundary conditions
+domain.set_boundary({'left': Br, 
+                     'right': Br, 
+                     'top': Br, 
+                     'bottom': 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)  ) )
+
+#------------------------------------------------------------------------------
+#
+# Evolve system through time
+#
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep=0.3, finaltime=40.0):
+    print domain.timestepping_statistics()
+
Index: /trunk/anuga_work/development/gareth/tests/urban_flow/plotme.py
===================================================================
--- /trunk/anuga_work/development/gareth/tests/urban_flow/plotme.py	(revision 8383)
+++ /trunk/anuga_work/development/gareth/tests/urban_flow/plotme.py	(revision 8383)
@@ -0,0 +1,48 @@
+### Python code to do plotting of this experiment.
+
+import util
+import numpy
+from matplotlib import pyplot as pyplot
+
+# Define gauge points, following Soares Frazao et al (2008)
+offset=numpy.array([7.7,+1.8])
+G1=numpy.array([2.65,1.15]) + offset
+G2=numpy.array([2.65, -0.60]) + offset
+G3=numpy.array([4.00, 1.15]) + offset
+G4=numpy.array([4.00,-0.8]) + offset
+G5=numpy.array([5.20,0.3]) + offset
+G6=numpy.array([-1.87,1.10]) + offset
+
+p=util.get_output('urban_flow0p1.sww')
+
+G1ind = ((p.x-G1[0])**2 + (p.y-G1[1])**2).argmin()
+pyplot.plot(p.time, p.stage[:,G1ind], 'o-')
+pyplot.savefig('G1.png')
+
+pyplot.clf()
+G2ind = ((p.x-G2[0])**2 + (p.y-G2[1])**2).argmin()
+pyplot.plot(p.time, p.stage[:,G2ind], 'o-')
+pyplot.savefig('G2.png')
+
+
+pyplot.clf()
+G3ind = ((p.x-G3[0])**2 + (p.y-G3[1])**2).argmin()
+pyplot.plot(p.time, p.stage[:,G3ind], 'o-')
+pyplot.savefig('G3.png')
+
+
+pyplot.clf()
+G4ind = ((p.x-G4[0])**2 + (p.y-G4[1])**2).argmin()
+pyplot.plot(p.time, p.stage[:,G4ind], 'o-')
+pyplot.savefig('G4.png')
+
+
+pyplot.clf()
+G5ind = ((p.x-G5[0])**2 + (p.y-G5[1])**2).argmin()
+pyplot.plot(p.time, p.stage[:,G5ind], 'o-')
+pyplot.savefig('G5.png')
+
+pyplot.clf()
+G6ind = ((p.x-G6[0])**2 + (p.y-G6[1])**2).argmin()
+pyplot.plot(p.time, p.stage[:,G6ind], 'o-')
+pyplot.savefig('G6.png')
