Index: /trunk/anuga_core/source/anuga/utilities/plot_utils.py
===================================================================
--- /trunk/anuga_core/source/anuga/utilities/plot_utils.py	(revision 9353)
+++ /trunk/anuga_core/source/anuga/utilities/plot_utils.py	(revision 9354)
@@ -1039,5 +1039,5 @@
 
 
-            if(bounding_polygon is not None):
+            if ( (bounding_polygon is not None) and (len(cut_points)>0)):
                 # Cut the points outside the bounding polygon
                 gridq[cut_points]= numpy.nan
@@ -1060,22 +1060,48 @@
     return
 
-def plot_triangles(p, adjustLowerLeft=False):
+def plot_triangles(p, adjustLowerLeft=False, values=None, color='k'):
     """ Add mesh triangles to a pyplot plot
-    """
+        
+       @param p = object holding sww vertex information (from util.get_output)
+       @param adjustLowerLeft = if TRUE, use spatial coordinates, otherwise use ANUGA internal coordinates     
+       @param values = list or array of length(p.x), or None. All triangles are assigned this value (for face plotting colors).
+       @param color = edge color for polygons (using matplotlib.colors notation)
+    """
+    import matplotlib
     from matplotlib import pyplot as pyplot
+    from matplotlib.collections import PolyCollection
     #
     x0=p.xllcorner
     x1=p.yllcorner 
-    #
+
+    # Get current figure information
+    #ax = pyplot.axes()
+
+    ##
+    vertices = []
     for i in range(len(p.vols)):
         k1=p.vols[i][0]
         k2=p.vols[i][1]
         k3=p.vols[i][2]
-        if(not adjustLowerLeft):
-            pyplot.plot([p.x[k1], p.x[k2], p.x[k3], p.x[k1]], [p.y[k1], p.y[k2], p.y[k3], p.y[k1]],'-',color='black')
+        if not adjustLowerLeft:
+            vertices.append([ [p.x[k1], p.y[k1]], [p.x[k2], p.y[k2]], [p.x[k3], p.y[k3]] ])
         else:
-            pyplot.plot([p.x[k1]+x0, p.x[k2]+x0, p.x[k3]+x0, p.x[k1]+x0], [p.y[k1]+x1, p.y[k2]+x1, p.y[k3]+x1, p.y[k1]+x1],'-',color='black')
-        #pyplot.plot([p.x[k3], p.x[k2]], [p.y[k3], p.y[k2]],'-',color='black')
-        #pyplot.plot([p.x[k3], p.x[k1]], [p.y[k3], p.y[k1]],'-',color='black')
+            vertices.append([ [p.x[k1]+x0, p.y[k1]+y0], [p.x[k2]+x0, p.y[k2]+y0], [p.x[k3]+x0, p.y[k3]+y0] ])
+      
+    if values is None: 
+        all_poly = PolyCollection( vertices, array = numpy.zeros(len(p.vols)), 
+            edgecolors=color)
+        all_poly.set_facecolor('none')
+    else:
+        assert len(values)==len(p.vols), 'len(values) must either be 1, or the same as len(p.vols)'
+        all_poly = PolyCollection( vertices, array = values, cmap = matplotlib.cm.jet, edgecolors=color)
+
+    pyplot.gca().add_collection(all_poly)
+    #    if(not adjustLowerLeft):
+    #        pyplot.plot([p.x[k1], p.x[k2], p.x[k3], p.x[k1]], [p.y[k1], p.y[k2], p.y[k3], p.y[k1]],'-',color=color)
+    #    else:
+    #        pyplot.plot([p.x[k1]+x0, p.x[k2]+x0, p.x[k3]+x0, p.x[k1]+x0], [p.y[k1]+x1, p.y[k2]+x1, p.y[k3]+x1, p.y[k1]+x1],'-',color=color)
+    #    #pyplot.plot([p.x[k3], p.x[k2]], [p.y[k3], p.y[k2]],'-',color='black')
+    #    #pyplot.plot([p.x[k3], p.x[k1]], [p.y[k3], p.y[k1]],'-',color='black')
 
 def find_neighbours(p,ind):
Index: /trunk/anuga_core/validation_tests/case_studies/merewether/Observations/ObservationPoints.csv
===================================================================
--- /trunk/anuga_core/validation_tests/case_studies/merewether/Observations/ObservationPoints.csv	(revision 9354)
+++ /trunk/anuga_core/validation_tests/case_studies/merewether/Observations/ObservationPoints.csv	(revision 9354)
@@ -0,0 +1,6 @@
+x, y, ID, stage (Draft Report), stage (ARR Report Final), TUFLOW (ARR Report Final)
+382373.514531995810103, 6354387.837079666554928, 4, 23.0, 23.01, 22.77
+382354.610297574603464, 6354365.208338395692408, 3, 23.1, 23.14, 23.11
+382424.399931652704254, 6354478.333491845987737, 0, 20.0, 19.98, 20.08
+382509.713526756153442, 6354548.220816057175398, 1, 18.4, 18.38, 18.36
+382339.416016335249878, 6354297.836651652120054, 2, 23.5, 23.36, 23.56
Index: /trunk/anuga_core/validation_tests/case_studies/merewether/Road/RoadPolygon.csv
===================================================================
--- /trunk/anuga_core/validation_tests/case_studies/merewether/Road/RoadPolygon.csv	(revision 9354)
+++ /trunk/anuga_core/validation_tests/case_studies/merewether/Road/RoadPolygon.csv	(revision 9354)
@@ -0,0 +1,96 @@
+382213.817011042,6354235.24761079
+382229.840596069,6354242.9389316
+382252.487262905,6354257.46698202
+382274.492986341,6354272.63597585
+382291.584810369,6354288.44591307
+382315.513364008,6354314.29729692
+382330.896005633,6354327.54346054
+382343.501225854,6354339.93503296
+382357.388332877,6354355.74497018
+382377.89852171,6354376.25515902
+382386.871729325,6354384.80107103
+382409.091100561,6354391.21050504
+382435.583427804,6354398.47453025
+382447.120409023,6354402.74748626
+382442.847453016,6354420.69390149
+382436.651666806,6354438.21302112
+382429.387641594,6354465.13264396
+382401.827075349,6354457.01402755
+382381.316886516,6354451.03188914
+382375.334748106,6354462.14157476
+382383.88066012,6354466.41453076
+382398.836006144,6354464.27805276
+382414.218647769,6354470.26019117
+382429.601289394,6354477.09692078
+382447.975000224,6354481.79717239
+382465.066824252,6354486.497424
+382479.594874676,6354489.9157888
+382477.458396672,6354501.45277002
+382472.330849464,6354513.41704684
+382469.76707586,6354526.23591486
+382469.339780259,6354532.21805327
+382480.449465877,6354535.20912248
+382499.677767909,6354539.90937408
+382511.428396928,6354542.90044329
+382529.161164357,6354546.53245589
+382543.048271379,6354551.4463553
+382540.911793376,6354562.12874532
+382533.861415964,6354589.68931156
+382525.315503951,6354622.59107282
+382521.683491345,6354633.05981504
+382500.53235911,6354626.65038102
+382464.853176452,6354617.46352561
+382479.808522476,6354633.70075844
+382519.547013341,6354640.53748805
+382516.769591937,6354651.43352586
+382527.879277555,6354651.43352586
+382529.588459958,6354644.16950065
+382551.594183393,6354649.93799126
+382566.335881617,6354651.64717367
+382583.641353446,6354651.86082147
+382583.427705645,6354649.08340006
+382529.802107758,6354637.11912324
+382541.339088977,6354600.58534938
+382549.24405759,6354567.46994033
+382553.303365796,6354554.86472011
+382603.296951078,6354568.32453153
+382606.928963684,6354561.70144972
+382556.935378402,6354547.3870471
+382582.294784283,6354447.5086605
+382572.235368431,6354444.48452806
+382545.949177293,6354544.09118676
+382537.707076371,6354541.61855649
+382504.378019516,6354531.57710987
+382480.876761478,6354526.44956266
+382478.953931275,6354520.25377645
+382483.654182883,6354504.44383922
+382487.499843289,6354489.0611976
+382483.226887282,6354483.50635479
+382469.126132459,6354477.52421638
+382453.316195233,6354474.31949938
+382437.933553608,6354467.69641757
+382441.138270613,6354456.58673195
+382447.761352424,6354430.0944047
+382454.598082035,6354407.87503347
+382494.3365729,6354418.77107129
+382525.742799551,6354429.6671091
+382543.048271379,6354433.94006511
+382573.81355463,6354443.12692053
+382605.219781281,6354450.81824134
+382607.142611484,6354436.29019091
+382553.730661397,6354423.68497069
+382515.915000735,6354411.07975047
+382483.867830683,6354398.68817805
+382458.871038042,6354391.63780064
+382427.678459191,6354382.87824083
+382409.067576944,6354377.92124514
+382395.844936939,6354371.12761181
+382378.753112911,6354359.59063059
+382347.560534061,6354324.97968693
+382322.77738922,6354302.1193723
+382306.753804194,6354284.81390047
+382289.021036765,6354266.65383744
+382265.306130926,6354250.20295681
+382231.763426272,6354230.76100698
+382218.730910451,6354223.28333397
+382213.817011042,6354235.24761079
Index: /trunk/anuga_core/validation_tests/case_studies/merewether/merewether.py
===================================================================
--- /trunk/anuga_core/validation_tests/case_studies/merewether/merewether.py	(revision 9353)
+++ /trunk/anuga_core/validation_tests/case_studies/merewether/merewether.py	(revision 9354)
@@ -7,4 +7,6 @@
 
 Water Research Laboratory, UNSW
+
+Edits by Steve Roberts & Gareth Davies, 2014
 """
 
@@ -22,12 +24,14 @@
 # Related major packages
 import anuga
-from anuga_parallel import myid, distribute, finalize
+from anuga_parallel import myid, distribute, finalize, barrier
 
 from anuga_parallel.parallel_operator_factory import Inlet_operator
+from anuga.utilities import quantity_setting_functions as qs
+from anuga import plot_utils as util
 
 # Application specific imports
 import project                 # Definition of file names and polygons
 
-verbose = project.verbose
+#verbose = project.verbose
 use_cache = project.use_cache
 
@@ -41,4 +45,8 @@
 meshname = 'merewether.msh'
 dem_name = 'topography1.dem'
+
+args = anuga.get_args()
+alg = args.alg
+verbose = args.verbose
 
 if myid == 0:
@@ -60,18 +68,44 @@
 holes_res = 1.0
 interior_regions = [[project.poly_merewether, merewether_res]]
-holes = project.holes
+
+# Either use houses as holes, or as part of the mesh (with breaklines)
+houses_as_holes = False
+if houses_as_holes:
+    holes = project.holes
+    breaklines = []
+else:
+    # Houses as elevation
+    house_height = 3.0
+    holes = []
+    breaklines = project.holes
+    house_addition_poly_fun_pairs = []
+    for i in range(len(breaklines)):
+        breaklines[i] = breaklines[i] + [breaklines[i][0]]
+        house_addition_poly_fun_pairs.append( 
+            [ breaklines[i], house_height])
+    house_addition_poly_fun_pairs.append(['All', 0.])
+
+#------------------------------------------------------------------------------
+# Make domain
+#------------------------------------------------------------------------------
 
 if myid == 0:
-    domain = anuga.create_domain_from_regions(project.bounding_polygon,
-   		boundary_tags={'bottom': [0],
-   			       'right': [1],
-                 	       'top': [2],
-                 	       'left': [3]},
-  				maximum_triangle_area=remainder_res,
-				interior_holes=holes,
-				mesh_filename=meshname,
-   				interior_regions=interior_regions,
-   				use_cache=use_cache,
-   				verbose=verbose)
+    mesh = anuga.create_mesh_from_regions(
+            project.bounding_polygon,
+            boundary_tags={'bottom': [0],
+                           'right': [1],
+                           'top': [2],
+                           'left': [3]},
+            maximum_triangle_area=remainder_res,
+            interior_holes=holes,
+            breaklines=breaklines,
+            filename=meshname,
+            interior_regions=interior_regions,
+            use_cache=use_cache,
+            verbose=verbose)
+
+    domain = anuga.create_domain_from_file(meshname)
+    
+    domain.set_flow_algorithm(alg)
 
     #------------------------------------------------------------------------------
@@ -79,9 +113,34 @@
     #------------------------------------------------------------------------------
     domain.set_quantity('stage', 0.0)
-    domain.set_quantity('friction', 0.02)
-    domain.set_quantity('elevation', filename='topography1.pts',
-			   	          use_cache=use_cache,
-	              			  verbose=verbose)
-
+
+    # Friction -- 2 options
+    variable_friction = True
+    if not variable_friction:
+        # Constant friction
+        domain.set_quantity('friction', 0.02)
+
+    else:
+        # Set friction to 0.02 on roads, 0.04 elsewhere
+        road_polygon = anuga.read_polygon('Road/RoadPolygon.csv')
+        friction_function = qs.composite_quantity_setting_function(
+            [ [road_polygon, 0.02], ['All', 0.04] ], 
+            domain)
+        domain.set_quantity('friction', friction_function)
+
+    # Elevation
+    if houses_as_holes:
+        domain.set_quantity('elevation', filename='topography1.pts',
+                              use_cache=use_cache,
+                                  verbose=verbose)
+
+    else:
+        domain.set_quantity('elevation', filename='topography1.pts',
+                              use_cache=use_cache,
+                              verbose=verbose, location='vertices')
+        # Add house_height inside houses    
+        house_addition_function = qs.composite_quantity_setting_function(
+            house_addition_poly_fun_pairs, domain)
+        domain.add_quantity('elevation', house_addition_function, 
+                            location='centroids')
 else:
     domain = None
@@ -93,4 +152,6 @@
 domain = distribute(domain)
 
+domain.set_store_vertices_uniquely()
+
 #------------------------------------------------------------------------------
 # Setup computational domain
@@ -98,5 +159,5 @@
 domain.set_name('merewether_1m') # Name of sww file
 domain.set_datadir('.') # Store sww output here
-domain.set_minimum_storable_height(0.001) # Store only depth > 1cm
+#domain.set_minimum_storable_height(0.001) # Store only depth > 1cm
 
 
@@ -112,5 +173,5 @@
 
 domain.set_boundary({'interior': Br,
-		     'bottom':   Br,
+                     'bottom':   Br,
                      'right':    Bt, # outflow
                      'top':      Bt, # outflow
@@ -121,7 +182,7 @@
 
 #fixed_inflow = anuga.Inflow(domain,
-#			center=(382300.0,6354290.0),
-#			radius=15.00,
-#			rate=19.7)
+#           center=(382300.0,6354290.0),
+#           radius=15.00,
+#           rate=19.7)
 #domain.forcing_terms.append(fixed_inflow)
 #hydrograph = anuga.Inflow(center=(382300.0,6354290.0),radius=30.0,rate=anuga.file_function('test_hydrograph2.tms', quantities=['hydrograph']) 
@@ -132,10 +193,24 @@
 #------------------------------------------------------------------------------
 
-for t in domain.evolve(yieldstep=10,finaltime=2000):
-	if myid == 0: print domain.timestepping_statistics()
+for t in domain.evolve(yieldstep=10,finaltime=1000):
+    if myid == 0: print domain.timestepping_statistics()
 
 
 domain.sww_merge()
 
+#------------------------------------------------------------------------------
+# Make geotif output
+#------------------------------------------------------------------------------
+
+barrier()
+if myid==0:
+    util.Make_Geotif('merewether_1m.sww', 
+                     output_quantities=['depth', 'velocity', 
+                                        'friction', 'elevation'],
+                     myTimeStep='last',
+                     CellSize=1.0,
+                     EPSG_CODE=32756,
+                     bounding_polygon=project.bounding_polygon,
+                     k_nearest_neighbours=1)
 
 finalize()
Index: /trunk/anuga_core/validation_tests/case_studies/merewether/plot_results.py
===================================================================
--- /trunk/anuga_core/validation_tests/case_studies/merewether/plot_results.py	(revision 9354)
+++ /trunk/anuga_core/validation_tests/case_studies/merewether/plot_results.py	(revision 9354)
@@ -0,0 +1,99 @@
+from anuga.utilities import plot_utils as util
+from matplotlib import pyplot as pyplot
+import numpy
+
+p=util.get_output('merewether_1m.sww')
+p2=util.get_centroids(p)
+# Time index at last time
+tindex = len(p2.time)-1
+
+# Get nearest wet points to 'point observations'
+point_observations = numpy.genfromtxt(
+    'Observations/ObservationPoints.csv',
+    delimiter=",",skip_header=1)
+
+nearest_points = point_observations[:,0]*0. - 1
+for i in range(len(nearest_points)):
+    # Compute distance of ANUGA points to observation, and
+    # if the ANUGA point is dry then add a large value
+    # Then find index of minimum
+    n = ( (p2.x+p2.xllcorner-point_observations[i,0])**2 + \
+          (p2.y+p2.yllcorner-point_observations[i,1])**2 + \
+          (p2.stage[tindex,:] <= p2.elev)*1.0e+06).argmin()
+    nearest_points[i] = n
+
+f = open('Stage_point_comparison.csv','w')
+f.writelines( 'Field Observation, ANUGA, TUFLOW, ANUGA minus Field Observation, ANUGA minus TUFLOW \n' )
+for i in range(len(nearest_points)):
+    po = point_observations[i,-2]
+    tu = point_observations[i,-1]
+    anuga_data = p2.stage[tindex, nearest_points.tolist()[i]]
+    newline = str(round(po,2)) + ', ' + str(round(anuga_data,2)) + ', ' + str(tu) + ', ' + \
+          str(round(anuga_data - po,2)) + ', ' + str(round(anuga_data - tu,2)) + '\n'
+    f.writelines(newline)
+
+f.flush()
+f.close()
+
+## Plot transect 1 [need to guess appropriate end points as these are not so
+## clear from the report]
+xx=util.near_transect(p2,[103, 100.], [130.,80.],tol=0.5)
+xx2=xx[0]
+
+pyplot.clf()
+pyplot.figure(figsize=(16,10.5))
+pyplot.subplot(121)
+pyplot.scatter(p2.x, p2.y, c=p2.elev,edgecolors='none')
+# Add nice elevation data
+colVals = numpy.maximum(numpy.minimum(p2.elev, 25.), 19.)
+util.plot_triangles(p, values = colVals, color='none')
+
+pyplot.gca().set_aspect('equal')
+pyplot.scatter(p2.x[xx2],p2.y[xx2],color='green')
+pyplot.xlim( (40., 160.))
+pyplot.ylim( (0.,140.))
+pyplot.title('Transect points in green')
+
+pyplot.subplot(222)
+pyplot.scatter(xx[1],p2.vel[tindex,xx[0]])
+#pyplot.xlim(0,25)
+pyplot.title('Final flow speed along the transect')
+
+pyplot.subplot(224)
+pyplot.scatter(xx[1],p2.stage[tindex,xx[0]]-p2.elev[xx[0]])
+#pyplot.xlim(0,25)
+pyplot.title('Final depth along the transect')
+pyplot.savefig('Transect1.png', bbox_inches='tight')
+
+pyplot.clf()
+
+# Velocity vector plot
+pyplot.figure(figsize=(16,22))
+pyplot.scatter(p2.x,p2.y,c=(p2.elev>24.),edgecolors='none', s=0.2)
+pyplot.gca().set_aspect('equal')
+pyplot.xlim((100,180))
+pyplot.ylim((100,210))
+#k=range(0,len(p2.x),2) # Thin out the vectors for easier viewing
+colVals = numpy.maximum(numpy.minimum(p2.elev, 25.), 19.)
+util.plot_triangles(p, values = colVals, color='white')
+k = range(len(p2.x))
+# Thin out the triangles
+#k = (((10.*(p2.x - p2.x.round())).round()%2 == 0.0)*((10.*(p2.y - p2.y.round())).round()%2 == 0.0)).nonzero()[0]
+pyplot.quiver(p2.x[k],p2.y[k],p2.xvel[tindex,k], p2.yvel[tindex,k],
+              scale_units='xy',units='xy',width=0.1,
+              color='black',scale=1.0)
+pyplot.savefig('velocity_stationary.png',dpi=100, bbox_inches='tight')
+
+
+pyplot.clf()
+pyplot.figure(figsize=(16,22))
+froude_number = p2.vel[tindex]/(numpy.maximum(p2.height[tindex], 1.0e-03)*9.8)**0.5
+froude_category = (froude_number>1.).astype(float) + (froude_number > 0.).astype(float)
+pyplot.scatter(p2.x,p2.y,edgecolors='none', s=0.2)
+pyplot.gca().set_aspect('equal')
+util.plot_triangles(p, values = froude_category, color='none')
+pyplot.xlim((p.x.min(), p.x.max()))
+pyplot.ylim((p.y.min(), p.y.max()))
+pyplot.title("Froude Number zones: 0, (0,1], or >1")
+pyplot.savefig('froudeNumber.png',dpi=100,bbox_inches='tight')
+
Index: /trunk/anuga_core/validation_tests/case_studies/merewether/produce_results.py
===================================================================
--- /trunk/anuga_core/validation_tests/case_studies/merewether/produce_results.py	(revision 9354)
+++ /trunk/anuga_core/validation_tests/case_studies/merewether/produce_results.py	(revision 9354)
@@ -0,0 +1,5 @@
+import os
+
+os.system('mpirun -np 6 python merewether.py -alg DE0')
+os.system('python plot_results.py')
+os.sytemt('pdflatex report.tex')
Index: /trunk/anuga_core/validation_tests/case_studies/merewether/project.py
===================================================================
--- /trunk/anuga_core/validation_tests/case_studies/merewether/project.py	(revision 9353)
+++ /trunk/anuga_core/validation_tests/case_studies/merewether/project.py	(revision 9354)
@@ -19,5 +19,5 @@
 # bounding polygon for study area
 bounding_polygon = anuga.read_polygon('extent.csv')
-print 'Area of bounding polygon', anuga.polygon_area(bounding_polygon)/1000000.0
+print 'Area of bounding polygon (km^2)', anuga.polygon_area(bounding_polygon)/1.0e+06
 
 
Index: /trunk/anuga_core/validation_tests/case_studies/merewether/results.tex
===================================================================
--- /trunk/anuga_core/validation_tests/case_studies/merewether/results.tex	(revision 9354)
+++ /trunk/anuga_core/validation_tests/case_studies/merewether/results.tex	(revision 9354)
@@ -0,0 +1,78 @@
+\documentclass{article}
+\usepackage{datatool}
+\usepackage{graphicx}
+\usepackage{datatool}
+
+\begin{document}
+
+\title{ANUGA applied to the ARR Project 15 Urban Flood Modelling test case}
+
+\maketitle
+
+\section{Introduction}
+This is a model of the Pasha Bulker flood (June 2007) in part of Merewether,
+Newcastle, Australia. It was developed by the UNSW Water Reseach Laboratory (WRL),
+and various information on the event is reported in `ARR Project 15: Two
+dimensional simulations in urban areas, Representation of buildings in 2D
+numerical flood models'. This chapter should be consulted for details on this
+test case, and reports on the performance of other models. 
+
+The ANUGA model was initially developed by WRL (who gave us permission to
+include it here) and further refined by ANUGA developers to include variable
+friction (with roads having a manning's n of 0.02, and everywhere else having
+0.04), an inlet operator, and options for either using buildings as topography, or
+buildings cut out of the mesh.
+
+\section{Results}
+Table~\ref{tab:comparison} compares the observed peak stage during the flood to
+ANUGA, and also values for a TUFLOW model which are reported in the ARR Project
+15 report. The ANUGA values should be similar to the field observations, and
+even more-so to the TUFLOW numerical model values. This suggests both models
+are giving a similar solution, with the main source of errors being related to
+the input data and the shallow water assumptions, rather than the numerical
+scheme.
+
+\DTLloaddb{comparison}{Stage_point_comparison.csv}
+\begin{table}
+\caption{Comparison of peak stage field observations, the ANUGA model, and a
+TUFLOW model (developed by WRL). See the ARR project 15 report for more
+information}
+\label{tab:comparison}
+\DTLdisplaydb{comparison}
+\end{table}
+
+Figure~\ref{fig:stationary_vel} shows a vector plot of the velocities around some buildings. It should
+show a similar pattern to Figure 30 in the ARR Project 15 report.
+\begin{figure}
+\includegraphics[width=\textwidth]{velocity_stationary.png}
+\caption{Velocity vectors around some buildngs at the end of the simulation
+(where it should have reached stationary state). Compare with Figure 30 in the ARR Project 15 report (the results should be similar)}
+\label{fig:stationary_vel}
+\end{figure}
+
+Figure~\ref{fig:transect1} shows a transect of velocity and depth between 2 buildings (Transect
+1 in the ARR Report).  It should show a similar pattern to the TUFLOW and
+MIKE21 models described in the ARR Report (their Figures 27 and 28).  In this
+site the numerical models show some differences to the physical model data,
+which may reflect the limitations both of the physical model, and of the
+shallow water assumptions (see ARR Report for further discussion).
+\begin{figure}
+\includegraphics[width=\textwidth]{Transect1.png}
+\caption{Velocity and depth along Transect 1. Compare with Figures 27 and 28 
+in the ARR Project 15 report (the results should be similar to TUFLOW and
+MIKE21)}
+\label{fig:transect1}
+\end{figure}
+
+Figure~\ref{fig:froude} shows the froude number (as either 0, $(0-1]$ or $>1$).  It
+should show a similar pattern to Figure 29 in the ARR report, although the
+colors are different. (FIXME: At the moment the zones can only be determined by
+looking at the ARR Report -- the plot needs a legend etc).
+\begin{figure}
+\center
+\includegraphics[width=0.7\textwidth]{froudeNumber.png}
+\caption{Froude number (0, $(0-1]$, or $>1$). The largest category area is zero, the smallest category area is $>1$. Compare with Figure 29 in ARR report (FIXME: FIX THE PLOT)}
+\label{fig:froude}
+\end{figure}
+
+\end{document}
