Index: anuga_work/development/dam_test_from_brad_2007/dam.csv
===================================================================
--- anuga_work/development/dam_test_from_brad_2007/dam.csv	(revision 4219)
+++ anuga_work/development/dam_test_from_brad_2007/dam.csv	(revision 4219)
@@ -0,0 +1,4 @@
+0,0
+6,0
+6,10
+0,10
Index: anuga_work/development/dam_test_from_brad_2007/dam_sample.py
===================================================================
--- anuga_work/development/dam_test_from_brad_2007/dam_sample.py	(revision 4219)
+++ anuga_work/development/dam_test_from_brad_2007/dam_sample.py	(revision 4219)
@@ -0,0 +1,105 @@
+"""
+Dam break!
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import os
+import time
+from shutil import copy
+from os import mkdir, access, F_OK, sep, environ, getenv, getcwd
+from os.path import expanduser
+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.utilities.polygon import read_polygon, plot_polygons, \
+                                    polygon_area, is_inside_polygon
+from anuga.utilities.polygon import Polygon_function
+from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 
+from anuga.geospatial_data.geospatial_data import *
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+length = 20.
+width = 10.
+dx = dy = 2           # Resolution: Length of subdivisions on both axes
+
+points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
+                                               len1=length, len2=width)
+
+domain = Domain(points, vertices, boundary)   # Create domain
+domain.set_name('dam_sample')                 # Output name
+domain.set_store_vertices_uniquely(True)      # Look at the 'real' triangles
+domain.beta_h = 0
+#domain.alpha_balance = 100000000
+
+
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+# Polygons
+study = [[0,0], [20,0], [20,10], [0,10], [0,0]]
+#dam = [[0,0], [6,0], [6,10], [0,10], [0,0]]
+poly_dam = read_polygon('dam.csv')
+#print poly_dam
+
+def dam1(x,y):
+    z = -x/10                 #initialise with slope (this is required to create the vector z?)
+    X = len(x)                  #find the number of vector elements
+    for i in range(X):
+        if 6<= x[i]<=7 or x[i]<1:    #dam wall
+                z[i] = 2
+        if 6<= x[i] <=7 and 4<= y[i] <=6: #remove unwanted wall
+                z[i] = -x[i]/10
+    return z
+
+def water(x,y):
+    z = -x/10                 #initialise with slope (this is required to create the vector z)
+    X = len(x)                #find the number of vector elements
+
+    for i in range(X):
+        if x[i]<=6:    #dam
+            z[i]= 2
+            
+    return z
+"""
+# maunal definition
+
+domain.set_quantity('elevation', water)                # Use function for elevation
+domain.set_quantity('friction', 0.01)                       # Constant friction 
+domain.set_quantity('stage', expression='elevation')        # take elevation
+domain.set_quantity('elevation', dam1)                # Use function for elevation
+"""
+
+# polygon definition
+domain.set_quantity('elevation', dam1)                # Use function for elevation
+domain.set_quantity('friction', 0.01)                       # Constant friction 
+domain.set_quantity('stage', Polygon_function( [(study, -1),(poly_dam, 2)] )) 
+
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([1.5, 0, 0])         # Inflow
+Bo = Dirichlet_boundary([-20/10, 0, 0])         # Outflow
+Br = Reflective_boundary(domain)             # Solid reflective wall
+
+#domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 0.2, finaltime = 10.0):
+    domain.write_time()
+
