Index: /trunk/anuga_core/demos/operators/run_bed_shear_erosion.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_bed_shear_erosion.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_bed_shear_erosion.py	(revision 9507)
@@ -0,0 +1,126 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+
+import numpy as num
+
+#===============================================================================
+# Setup Functions
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+def topography(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+    print ' Create topography....'
+    z = -x/100
+
+    N = len(x)
+    for i in range(N):
+        # Step
+        if 2 < x[i] < 4:
+            z[i] += 0.4 - 0.05*y[i]
+
+        # Permanent pole
+        if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
+            z[i] += 1
+
+        # Dam
+        #if 12 < x[i] < 13:
+        #    z[i] += 1.4
+    # Sloping Embankment Across Channel
+        if 10.6 < x[i] < 12.1:
+            # Cut Out Segment for Culvert face                
+            z[i] +=  1.0*(x[i] -10.6)    # Sloping Segment  U/S Face
+        if 12.0 < x[i] < 13.0:
+           z[i] +=  1.4                    # Flat Crest of Embankment
+        if 12.9 < x[i] < 13.7:
+            z[i] +=  1.4-2.0*(x[i] -12.9) # Sloping D/S Face
+            
+    return z
+
+
+
+#===============================================================================
+# Setup and Run Model
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+print ' Set up Domain first...'
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_flow_algorithm('DE0')
+domain.set_name() # Output name
+domain.set_store_vertices_uniquely(True)
+
+print domain.statistics()
+
+domain.set_quantities_to_be_stored({'elevation': 2,
+                                    'stage': 2,
+                                    'xmomentum': 2,
+                                    'ymomentum': 2})
+
+
+
+domain.set_quantity('elevation', topography)           # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([1.5, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Setup erosion operator in the middle of dam
+#------------------------------------------------------------------------------
+print 'Set up Erosion Area to test...'
+
+from anuga import Bed_shear_erosion_operator
+polygon1 = [ [10.6, 1.0], [13.7, 1.0], [13.7, 4.0], [10.6, 4.0] ]
+
+
+# create operator
+op1 = Bed_shear_erosion_operator(domain, threshold=2.0, shear_factor=1000.0, polygon=polygon1)
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep=0.1, finaltime=15.0):
+    domain.print_timestepping_statistics()
+    #domain.print_operator_timestepping_statistics()
+
+
+
+
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_change_elevation.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_change_elevation.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_change_elevation.py	(revision 9507)
@@ -0,0 +1,119 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import numpy
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_name()
+print domain.statistics()
+domain.set_quantities_to_be_stored({'elevation': 2,
+                                    'stage': 2,
+                                    'xmomentum': 2,
+                                    'ymomentum': 2})
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+def topography_dam(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+
+    z = -x/100
+
+    # Step
+    id = (2 < x) & (x < 4)
+    z[id] += 0.4 - 0.05*y[id]
+
+
+    # Permanent pole
+    #id = (x - 8)**2 + (y - 2)**2 < 0.4**2
+    #z[id] += 1
+
+
+    # Dam
+    id = (12 < x) & (x  < 13)
+    z[id] += 0.4
+
+
+            
+    return z
+
+def topography_dam_break(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+
+    z = -x/100
+    
+    # Step
+    id = (2 < x) & (x < 4)
+    z[id] += 0.4 - 0.05*y[id]
+
+
+    # Permanent pole
+    #id = (x - 8)**2 + (y - 2)**2 < 0.4**2
+    #z[id] += 1
+
+    # Dam with hole
+    id = (12 < x) & (x < 13) and (y > 3.0)
+    z[id] += 0.4
+
+    id = (12 < x) & (x < 13) and (y < 2.0)
+    z[id] += 0.4
+
+
+    return z
+
+
+domain.set_quantity('elevation', topography_dam)       # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+from anuga.operators.set_elevation import Set_elevation
+op1 = Set_elevation(domain, elevation=lambda x,y : -x/100, radius=1.0, center = (12.5,3.0))
+
+
+#dam_break = False
+
+for t in domain.evolve(yieldstep=0.1, finaltime=30.0):
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+    if numpy.allclose(t, 10.0):
+        print 'changing elevation'
+        op1()
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_collect_max_stage_operator.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_collect_max_stage_operator.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_collect_max_stage_operator.py	(revision 9507)
@@ -0,0 +1,95 @@
+"""Simple water flow example using ANUGA
+
+Water driven up a linear slope and time varying boundary,
+similar to a beach environment
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+import anuga
+
+
+from math import cos
+from numpy import zeros, float, where
+import numpy
+from time import localtime, strftime, gmtime
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1000.
+dy = dx
+L = 100000.
+W = 10*dx
+#W = dx
+
+# structured mesh
+domain = anuga.rectangular_cross_domain(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+domain.set_name() # based on script name
+
+print domain.mesh.statistics(nbins=50)
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+domain.set_flow_algorithm('2_0')
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+h0 = 10.0
+h1 = 0.0
+
+def elevation(x,y):
+
+    return where(x<=95000.0, 0.0, h0)
+
+def height(x,y):
+
+    return where(x<=50000.0, h0, h1)
+
+domain.set_quantity('elevation',elevation)
+domain.set_quantity('friction', 0.0)
+domain.add_quantity('stage', height)
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+from math import sin, pi, exp
+Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
+Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary 
+Bd = anuga.Dirichlet_boundary([1,0.,0.])    # Constant boundary values
+
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+#-----------------------------------------------------------------------------
+# Setup operators that will be applied each inner timestep
+#------------------------------------------------------------------------------
+from anuga.operators.collect_max_stage_operator import Collect_max_stage_operator
+max_operator = Collect_max_stage_operator(domain)
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 100.0, finaltime = 60*60.):
+    #print domain.timestepping_statistics(track_speeds=True)
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+
+# save the max_stage centroid data to a text file
+max_operator.save_centroid_data_to_csv()
+
+# Let's have a look at the max_stage
+#max_operator.plot_quantity()
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_flat_fill_slice_erosion.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_flat_fill_slice_erosion.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_flat_fill_slice_erosion.py	(revision 9507)
@@ -0,0 +1,134 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+
+import numpy as num
+
+#===============================================================================
+# Setup Functions
+#===============================================================================
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+def topography(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+    print ' Create topography....'
+    z = -x/100
+
+    N = len(x)
+    for i in range(N):
+        # Step
+        if 2 < x[i] < 4:
+            z[i] += 0.4 - 0.05*y[i]
+
+        # Permanent pole
+        if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
+            z[i] += 1
+
+        # Dam
+        #if 12 < x[i] < 13:
+        #    z[i] += 1.4
+    # Sloping Embankment Across Channel
+        if 10.6 < x[i] < 12.1:
+            # Cut Out Segment for Culvert face                
+            z[i] +=  1.0*(x[i] -10.6)    # Sloping Segment  U/S Face
+        if 12.0 < x[i] < 13.0:
+           z[i] +=  1.4                    # Flat Crest of Embankment
+        if 12.9 < x[i] < 13.7:
+            z[i] +=  1.4-2.0*(x[i] -12.9) # Sloping D/S Face
+            
+    return z
+
+
+
+#===============================================================================
+# Setup and Run Model
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+print ' Set up Domain first...'
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_flow_algorithm('DE1')
+domain.set_name('flat_fill_slice_erosion') # Output name
+print domain.statistics()
+
+domain.set_quantities_to_be_stored({'elevation': 2,
+                                    'stage': 2,
+                                    'xmomentum': 2,
+                                    'ymomentum': 2})
+
+
+
+domain.set_quantity('elevation', topography)           # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([1.5, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Setup erosion operator in the middle of dam
+#------------------------------------------------------------------------------
+print 'Set up Erosion Area to test...'
+
+from anuga import Flat_fill_slice_erosion_operator as Default_erosion_operator
+from anuga import Set_elevation_operator as Default_erosion_operator
+
+polygon1 = [ [10.6, 1.0], [13.7, 1.0], [13.7, 4.0], [10.6, 4.0] ]
+
+# setup interpolated function
+time = num.array([ 5.0 , 8.0])
+elev = num.array([ 1.4,  2.0])
+
+from scipy.interpolate import interp1d
+flevel = interp1d(time,elev)
+
+# create operator
+op1 = Default_erosion_operator(domain, elevation=flevel, polygon=polygon1)
+
+
+
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep=0.1, finaltime=15.0):
+    domain.print_timestepping_statistics()
+    #domain.print_operator_timestepping_statistics()
+
+
+
+
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_flat_slice_erosion.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_flat_slice_erosion.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_flat_slice_erosion.py	(revision 9507)
@@ -0,0 +1,133 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+
+import numpy as num
+
+#===============================================================================
+# Setup Functions
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+def topography(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+    print ' Create topography....'
+    z = -x/100
+
+    N = len(x)
+    for i in range(N):
+        # Step
+        if 2 < x[i] < 4:
+            z[i] += 0.4 - 0.05*y[i]
+
+        # Permanent pole
+        if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
+            z[i] += 1
+
+        # Dam
+        #if 12 < x[i] < 13:
+        #    z[i] += 1.4
+    # Sloping Embankment Across Channel
+        if 10.6 < x[i] < 12.1:
+            # Cut Out Segment for Culvert face                
+            z[i] +=  1.0*(x[i] -10.6)    # Sloping Segment  U/S Face
+        if 12.0 < x[i] < 13.0:
+           z[i] +=  1.4                    # Flat Crest of Embankment
+        if 12.9 < x[i] < 13.7:
+            z[i] +=  1.4-2.0*(x[i] -12.9) # Sloping D/S Face
+            
+    return z
+
+
+
+#===============================================================================
+# Setup and Run Model
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+print ' Set up Domain first...'
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_flow_algorithm('DE1')
+domain.set_name() # Output name
+print domain.statistics()
+
+domain.set_quantities_to_be_stored({'elevation': 2,
+                                    'stage': 2,
+                                    'xmomentum': 2,
+                                    'ymomentum': 2})
+
+
+
+domain.set_quantity('elevation', topography)           # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([1.5, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Setup erosion operator in the middle of dam
+#------------------------------------------------------------------------------
+print 'Set up Erosion Area to test...'
+
+from anuga import Flat_slice_erosion_operator
+polygon1 = [ [10.6, 1.0], [13.7, 1.0], [13.7, 4.0], [10.6, 4.0] ]
+
+# setup interpolated function
+time = num.array([ 2.0 , 10.0])
+elev = num.array([ 1.4,  -0.1])
+
+from scipy.interpolate import interp1d
+flevel = interp1d(time,elev)
+
+# create operator
+op1 = Flat_slice_erosion_operator(domain, elevation=flevel, polygon=polygon1)
+
+
+
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep=0.1, finaltime=15.0):
+    domain.print_timestepping_statistics()
+    #domain.print_operator_timestepping_statistics()
+
+
+
+
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_polygon_erosion.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_polygon_erosion.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_polygon_erosion.py	(revision 9507)
@@ -0,0 +1,176 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+
+
+
+#===============================================================================
+# Stick all the class and function definitions here
+#===============================================================================
+
+#-------------------------------------------------------------------------------
+# Here are the def's for defining quantities
+#-------------------------------------------------------------------------------
+def topography(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+
+    z = -x/100
+
+    # Step
+    id = (2 < x) & (x < 4)
+    z[id] += 0.4 - 0.05*y[id]
+
+    # Permanent pole
+    id = (x - 8)**2 + (y - 2)**2 < 0.4**2
+    z[id] += 1
+
+    # Dam
+    id = (12 < x) & (x < 13)
+    z[id] += 0.4
+
+    return z
+
+
+#-------------------------------------------------------------------------------
+# Inherit from erosion operator
+#-------------------------------------------------------------------------------
+from anuga.operators.erosion_operators import Polygonal_erosion_operator
+
+class My_polygon_erosion_operator(Polygonal_erosion_operator):
+    """
+    Local version of erosion confined to a polygon
+
+    """
+
+    def __init__(self, domain,
+                 threshold=0.0,
+                 base=0.0,
+                 polygon=None,
+                 verbose=False):
+
+
+        Polygonal_erosion_operator.__init__(self, domain, threshold, base, polygon, verbose)
+
+
+
+    def update_quantities(self):
+        """Update the vertex values of the quantities to model erosion
+        """
+        import numpy as num
+        
+        t = self.get_time()
+        dt = self.get_timestep()
+
+        updated = True
+
+        if self.indices is None:
+
+            #--------------------------------------
+            # Update all three vertices for each cell
+            #--------------------------------------
+            self.elev_v[:] = self.elev_v + 0.0
+
+        else:
+            ind = self.indices
+            m = num.sqrt(self.xmom_c[ind]**2 + self.ymom_c[ind]**2)
+            
+            if self.domain.get_using_discontinuous_elevation():
+                height = self.stage_c[ind] - self.elev_c[ind]
+                
+                m = num.where(m>self.threshold, m, 0.0)
+                self.elev_c[ind] = num.maximum(self.elev_c[ind] - m*dt, self.base)
+
+                self.stage_c[ind] = self.elev_c[ind] + height
+            else:
+            #--------------------------------------
+            # Update all three vertices for each cell
+            #--------------------------------------
+            
+            
+                m = num.vstack((m,m,m)).T
+                m = num.where(m>self.threshold, m, 0.0)
+                self.elev_v[ind] = num.maximum(self.elev_v[ind] - m*dt, self.base)
+                #num.maximum(self.elev_v[ind] - momentum*dt, Z)
+
+
+        return updated
+
+
+
+#===============================================================================
+# Now to the standard anuga model setup
+#===============================================================================
+
+
+#-------------------------------------------------------------------------------
+# Setup computational domain
+#-------------------------------------------------------------------------------
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_name() # Output name based on script
+domain.set_flow_algorithm('DE1')
+print domain.statistics()
+domain.set_store_vertices_uniquely(True)
+
+domain.set_quantities_to_be_stored({'elevation': 2,
+                                    'stage': 2,
+                                    'xmomentum': 2,
+                                    'ymomentum': 2})
+
+#-------------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+
+domain.set_quantity('elevation', topography)           # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#-------------------------------------------------------------------------------
+# Setup boundary conditions
+#-------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([0.5, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#-------------------------------------------------------------------------------
+# Setup erosion operator in the middle of dam
+#-------------------------------------------------------------------------------
+polygon1 = [ [12., 0.0], [13., 0.0], [13., 5.0], [12., 5.0] ]
+op1 = My_polygon_erosion_operator(domain, threshold=0.0, base=-0.1, polygon=polygon1)
+
+
+
+
+#-------------------------------------------------------------------------------
+# Evolve system through time
+#-------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep=0.2, finaltime=60.0):
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+
+
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_rate_operator.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_rate_operator.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_rate_operator.py	(revision 9507)
@@ -0,0 +1,146 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+import os
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_name() # Output name based on script name. You can add timestamp=True
+print domain.statistics()
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+def topography(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+
+    z = -x/100
+
+    # Step
+    id = (2 < x) & (x < 4)
+    z[id] += 0.4 - 0.05*y[id]
+
+
+    # Permanent pole
+    #id = (x - 8)**2 + (y - 2)**2 < 0.4**2
+    #z[id] += 1
+
+
+    # Dam
+    #id = (12 < x) & (x  < 13)
+    #z[id] += 0.4
+
+
+
+            
+    return z
+
+
+def pole_increment(x,y,t):
+    """This provides a small increment to a pole located mid stream
+    For use with variable elevation data
+    """
+
+    z = 0.0*x    
+
+    if t<10.0:
+        return z
+    
+
+    # Pole 1
+    id = (x - 12)**2 + (y - 3)**2 < 0.4**2
+    z[id] += 0.1
+
+
+    # Pole 2
+    id = (x - 14)**2 + (y - 2)**2 < 0.4**2
+    z[id] += 0.05
+
+
+    return z
+
+
+def pole(t):
+
+    if t<10:
+        return 0.0
+    elif t>15:
+        return 0.0
+    else:
+        return 0.05
+
+
+domain.set_quantity('elevation', topography)           # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+polygon1 = [ [10.0, 0.0], [11.0, 0.0], [11.0, 5.0], [10.0, 5.0] ]
+polygon2 = [ [12.0, 2.0], [13.0, 2.0], [13.0, 3.0], [12.0, 3.0] ]
+
+from anuga.operators.rate_operators import Rate_operator
+
+op1 = Rate_operator(domain, rate=lambda t: 10.0 if (t>=0.0) else 0.0, polygon=polygon2)
+op2 = Rate_operator(domain, rate=lambda t: 10.0 if (t>=0.0) else 0.0, radius=0.5, center=(10.0, 3.0))
+
+
+domain.set_starttime(-0.1)
+for t in domain.evolve(yieldstep=0.01, finaltime=0.0):
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+    stage = domain.get_quantity('stage')
+    elev  = domain.get_quantity('elevation')
+    height = stage - elev
+
+    print 'integral = ', height.get_integral()
+
+
+for t in domain.evolve(yieldstep=0.1, duration=5.0):
+
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+    stage = domain.get_quantity('stage')
+    elev  = domain.get_quantity('elevation')
+    height = stage - elev
+
+    print 'integral = ', height.get_integral()
+
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_rate_spatial_operator.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_rate_spatial_operator.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_rate_spatial_operator.py	(revision 9507)
@@ -0,0 +1,180 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import numpy
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+
+from anuga import indent
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_name() # Output name based on script
+print domain.statistics()
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+def topography(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+
+    z = -x/100
+
+    N = len(x)
+    for i in range(N):
+        # Step
+        if 2 < x[i] < 4:
+            z[i] += 0.4 - 0.05*y[i]
+
+        # Permanent pole
+        #if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
+        #    z[i] += 1
+
+
+#        # Dam
+#        if 12 < x[i] < 13 and y[i] > 3.0:
+#            z[i] += 0.4
+#
+#        if 12 < x[i] < 13 and y[i] < 2.0:
+#            z[i] += 0.4
+
+
+#        # Dam
+#        if 12 < x[i] < 13:
+#            z[i] += 0.4
+
+            
+    return z
+
+
+def pole_increment(x,y,t):
+    """This provides a small increment to a pole located mid stream
+    For use with variable elevation data
+    """
+
+    z = 0.0*x
+    
+    if t<10.0:
+        return z
+    
+    N = len(x)
+    for i in range(N):
+        # Pole 1
+        if (x[i] - 12)**2 + (y[i] - 3)**2 < 0.4**2:
+            z[i] += 0.1
+
+    for i in range(N):
+        # Pole 2
+        if (x[i] - 14)**2 + (y[i] - 2)**2 < 0.4**2:
+            z[i] += 0.05
+
+    return z
+
+
+def pole(t):
+
+    if t<10:
+        return 0.0
+    elif t>15:
+        return 0.0
+    else:
+        return 0.05
+
+
+domain.set_quantity('elevation', topography)           # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+polygon1 = [ [10.0, 0.0], [11.0, 0.0], [11.0, 5.0], [10.0, 5.0] ]
+polygon2 = [ [12.0, 2.0], [13.0, 2.0], [13.0, 3.0], [12.0, 3.0] ]
+
+
+from anuga.operators.rate_operators import Rate_operator
+
+op1 = Rate_operator(domain, rate=10.0, polygon=polygon2)
+
+area1 = numpy.sum(domain.areas[op1.indices])
+Q1 = 10.0*area1
+print 'op1 Q ',Q1
+
+op2 = Rate_operator(domain, rate=10.0, radius=0.5, center=(10.0, 3.0))
+
+area2 = numpy.sum(domain.areas[op2.indices])
+Q2 = 10.0*area2
+print 'op2 Q ',Q2
+
+
+def rain(x,y,t):
+    """Function to calculate "rain"
+    input x,y should be considered to be numpy arrays
+    abd t a scalar
+    """
+    if t<=4.0:
+        return (x+y)*1.0
+    else:
+        return 0.0*x
+
+factor = 1e-3
+op3 = Rate_operator(domain, rate = rain, factor=factor)
+Q3 = numpy.sum(op3.get_spatial_rate()*domain.areas)*factor
+
+
+#op3()
+#domain.fractional_step_operators.remove(op3)
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+accum = 0.0
+yieldstep = 0.1
+finaltime = 5.0
+for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+    stage = domain.get_quantity('stage')
+    elev  = domain.get_quantity('elevation')
+    height = stage - elev
+
+    print indent + 'Integral = ', height.get_integral()
+    print indent + 'Exact accumultion = ', accum
+    
+    dd = max(min(yieldstep,4.0-t),0.0)
+    accum += (Q1+Q2)*yieldstep + dd*Q3
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_set_depth_friction.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_set_depth_friction.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_set_depth_friction.py	(revision 9507)
@@ -0,0 +1,109 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+from anuga import Region
+from anuga import indent
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_name() # Output name
+print domain.statistics()
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+def topography(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+
+    z = -x/100
+
+    # Step
+    id = ( 2 < x ) & (x < 4)
+    z[id] +=  0.4 - 0.05*y[id]
+
+    # Permanent pole
+    id = (x - 8)**2 + (y - 2)**2 < 0.4**2
+    z[id] += 1
+
+    # Pole 2
+    #id =  (x - 14)**2 + (y - 3.5)**2 < 0.4**2
+    #z[id] += 1.0
+
+
+    return z
+
+
+
+domain.set_quantity('elevation', topography)           # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Setup operators which are applied each inner step
+#------------------------------------------------------------------------------
+from anuga.operators.set_friction_operators import Depth_friction_operator
+from scipy import interpolate
+
+h = [0.0,     0.2,   0.4,   0.8,   2.0, 99.0]
+n = [0.250, 0.060, 0.045, 0.035, 0.025, 0.025]
+
+friction = interpolate.interp1d(h, n, fill_value=0.025)
+op1 = Depth_friction_operator(domain,friction=friction)
+
+#p1 = [ [12.0, 2.5], [13.5, 2.5], [13.5, 4.0], [12.0, 4.0] ]
+#op2 = Depth_friction_operator(domain,
+#                                  friction = lambda h :10,
+#                                  polygon=p1)
+
+# Setup region for integrating quantities
+p2 = [ [8.0, 2.5], [9.5, 2.5], [9.5, 4.0], [8.0, 4.0] ]
+reg = Region(domain, polygon=p2)
+
+# Some useful aliases
+stage = domain.quantities['stage']
+elev = domain.quantities['elevation']
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep=0.1, finaltime=10.0):
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+    # let calculate the integral of height over a region
+    height = stage-elev
+    print indent+'Int_p2(h) = '+str(height.get_integral(region=reg))
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_set_elevation.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_set_elevation.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_set_elevation.py	(revision 9507)
@@ -0,0 +1,137 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+import numpy as num
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+
+
+#------------------------------------------------------------------------------
+# Stage and Elevation Functions
+#------------------------------------------------------------------------------
+def topography(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+
+
+    z = -x/100
+
+
+    # Step
+    id_s = (2 < x) & (x < 4)
+
+    z[id_s] += 0.4 - 0.05*y[id_s]
+
+    # Permanent pole
+    id_p = (x - 8)**2 + (y - 2)**2 < 0.4**2
+
+    z[id_p] += 1
+
+            
+    return z
+
+
+def pole_increment(x,y,t):
+    """This provides a small increment to a pole located mid stream
+    For use with variable elevation data
+    """
+
+
+    z = 0.0*x
+    
+
+    if t<10.0:
+        return z
+    
+
+    N = len(x)
+    for i in range(N):
+        # Pole 1
+        if (x[i] - 12)**2 + (y[i] - 3)**2 < 0.4**2:
+            z[i] += 0.1
+
+    for i in range(N):
+        # Pole 2
+        if (x[i] - 14)**2 + (y[i] - 2)**2 < 0.4**2:
+            z[i] += 0.05
+
+    return z
+
+
+def pole(x,y,t):
+	
+    z = topography(x,y)
+    if t<5:
+        return z
+    elif t>12:
+        return z
+    else:
+        return z+1.0
+
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_name() # Output name based on script name
+print domain.statistics()
+domain.set_quantities_to_be_stored({'elevation': 2,
+                                    'stage': 2})
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+domain.set_quantity('elevation', topography)           # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+from anuga import Set_elevation_operator
+op1 = Set_elevation_operator(domain, elevation=pole, radius=0.5, center = (12.0,3.0))
+
+from anuga import Set_elevation
+op2 = Set_elevation(domain, elevation = topography)
+
+growing = False
+shrinking = False
+
+done = False
+for t in domain.evolve(yieldstep=0.1, finaltime=20.0):
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+    if num.allclose(15.0, t):
+        op1.set_elevation(None)
+        op2()
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_set_elevation_polygon.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_set_elevation_polygon.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_set_elevation_polygon.py	(revision 9507)
@@ -0,0 +1,135 @@
+"""Simple water flow example using ANUGA
+
+Water flowing down a channel with a topography that varies with time
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+from anuga import rectangular_cross
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+length = 24.
+width = 5.
+dx = dy = 0.2 #.1           # 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)
+domain.set_name('set_elevation') # Output name
+print domain.statistics()
+domain.set_quantities_to_be_stored({'elevation': 2,
+                                    'stage': 2})
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+def topography_dam(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+
+    z = -x/100
+
+    N = len(x)
+    for i in range(N):
+        # Step
+        if 2 < x[i] < 4:
+            z[i] += 0.4 - 0.05*y[i]
+
+        # Permanent pole
+        if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
+            z[i] += 1
+
+
+        # Dam
+        if 12 < x[i] < 13:
+            z[i] += 0.4
+
+            
+    return z
+
+def topography_dam_break(x,y):
+    """Complex topography defined by a function of vectors x and y."""
+
+    z = -x/100
+
+    N = len(x)
+    for i in range(N):
+        # Step
+        if 2 < x[i] < 4:
+            z[i] += 0.4 - 0.05*y[i]
+
+        # Permanent pole
+        if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
+            z[i] += 1
+
+
+        # Dam
+        if 12 < x[i] < 13 and y[i] > 3.0:
+            z[i] += 0.4
+
+        if 12 < x[i] < 13 and y[i] < 2.0:
+            z[i] += 0.4
+
+    return z
+
+
+
+def pole(t):
+
+    if t<10:
+        return 0.0
+    elif t>15:
+        return 0.0
+    else:
+        return 0.05
+
+
+domain.set_quantity('elevation', topography_dam)       # elevation is a function
+domain.set_quantity('friction', 0.01)                  # Constant friction
+domain.set_quantity('stage', expression='elevation')   # Dry initial condition
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
+Br = Reflective_boundary(domain)              # Solid reflective wall
+Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+#from anuga.operators.set_elevation_operators import Circular_set_elevation_operator
+
+#op1 = Circular_set_elevation_operator(domain, elevation=pole, radius=0.5, center = (12.0,3.0))
+
+
+dam_break = False
+
+for t in domain.evolve(yieldstep=0.1, finaltime=40.0):
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+    if t >= 10 and not dam_break:
+        print 'changing elevation'
+
+        stage_c = domain.get_quantity('stage').centroid_values
+        elev_c =  domain.get_quantity('elevation').centroid_values
+        height_c = stage_c - elev_c
+        domain.set_quantity('elevation', topography_dam_break)
+        stage_c[:] = elev_c + height_c
+
+        dam_break = True
+
+
+
+
+
Index: /trunk/anuga_core/demos/operators/run_set_stage.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_set_stage.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_set_stage.py	(revision 9507)
@@ -0,0 +1,91 @@
+"""Simple test of setting stage in a circular region
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+import anuga
+
+
+from math import cos
+from numpy import zeros, float, where
+import numpy
+from time import localtime, strftime, gmtime
+
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 100.
+dy = dx
+L = 10000.
+W = L
+#W = dx
+
+# structured mesh
+domain = anuga.rectangular_cross_domain(int(L/dx), int(W/dy), L, W, (-L/2.0, -W/2.0))
+
+
+domain.set_name()                
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+domain.set_flow_algorithm(2.0)
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+h0 = 1000.0
+
+domain.set_quantity('elevation',0.0)
+domain.set_quantity('friction', 0.0)
+domain.add_quantity('stage', h0)
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Setup Operators
+#------------------------------------------------------------------------------
+from anuga.operators.set_stage_operator import Set_stage_operator
+
+import math
+
+stage1 = lambda t: h0 + 20.0 * math.sin(t/3.0)
+#cop1 = Set_stage_operator(domain, stage=stage1, center=(0.0, 0.0), radius=100.0 )
+
+stage2 = lambda t: h0 + 30.0# * math.sin(t/3.0)
+cop2 = Set_stage_operator(domain, stage=stage2, center=(2000.0, 1000.0), radius=100.0 )
+
+#print cop1.statistics()
+#print cop2.statistics()
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+## stage = domain.get_quantity('stage')
+## id = domain.get_triangle_containing_point([2000., 1150.])
+## times = []
+## stages = []
+
+for t in domain.evolve(yieldstep = 1.0, finaltime = 200.0):
+    #print domain.timestepping_statistics(track_speeds=True)
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+    #times.append(t)
+    #stages.append(stage.centroid_values[id])
+
+## import pylab as pl
+
+## pl.plot(times, stages)
+## pl.show()
Index: /trunk/anuga_core/demos/operators/run_set_w_uh_vh.py
===================================================================
--- /trunk/anuga_core/demos/operators/run_set_w_uh_vh.py	(revision 9507)
+++ /trunk/anuga_core/demos/operators/run_set_w_uh_vh.py	(revision 9507)
@@ -0,0 +1,79 @@
+"""Simple test of setting stage in a circular region
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+import anuga
+
+
+from math import cos
+from numpy import zeros, float, where
+import numpy
+from time import localtime, strftime, gmtime
+
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 100.
+dy = dx
+L = 10000.
+W = L
+#W = dx
+
+# structured mesh
+domain = anuga.rectangular_cross_domain(int(L/dx), int(W/dy), L, W, (-L/2.0, -W/2.0))
+
+
+domain.set_name()                
+
+#------------------------------------------------------------------------------
+# Setup Algorithm
+#------------------------------------------------------------------------------
+domain.set_flow_algorithm(2.0)
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+h0 = 1000.0
+
+domain.set_quantity('elevation',0.0)
+domain.set_quantity('friction', 0.0)
+domain.add_quantity('stage', h0)
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+#------------------------------------------------------------------------------
+# Setup Operators
+#------------------------------------------------------------------------------
+from anuga.operators.set_w_uh_vh_operator import Set_w_uh_vh_operator
+
+import math
+
+wuhvh1 = lambda t: [ h0 + 20.0 * math.sin(t/3.0) ,  20.0 * math.sin(t/3.0),  10.0 * math.sin(t/3.0)]
+cop1 = Set_w_uh_vh_operator(domain, w_uh_vh=wuhvh1, center=(0.0, 0.0), radius=100.0 )
+
+
+print cop1.statistics()
+#print cop2.statistics()
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 1.0, finaltime = 40.0):
+    #print domain.timestepping_statistics(track_speeds=True)
+    domain.print_timestepping_statistics()
+    domain.print_operator_timestepping_statistics()
+
+
Index: /trunk/anuga_core/demos/parallel/run_advection.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_advection.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_advection.py	(revision 9507)
@@ -0,0 +1,50 @@
+#import pdb
+#pdb.set_trace()
+
+#========================================================================
+from anuga.config import g, epsilon
+
+#from Numeric import allclose, array, zeros, ones, Float
+
+from anuga.advection import Domain, Transmissive_boundary, Dirichlet_boundary
+
+#from Numeric import array
+import numpy as num
+
+
+from anuga.interface import rectangular_cross
+
+#points, vertices, boundary = rectangular(60, 60)
+points, vertices, boundary = rectangular_cross(10, 10)
+
+#Create advection domain with direction (1,-1)
+domain = Domain(points, vertices, boundary, velocity=[1.0, 0.0])
+
+# Initial condition is zero by default
+
+#domain.initialise_visualiser(scale_z=1.0)
+#domain.visualise_range_z = 0.5
+#domain.visualise_color_stage = True
+
+
+#Boundaries
+T = Transmissive_boundary(domain)
+D = Dirichlet_boundary(num.array([1.0]))
+
+#turn on the visualisation
+rect = [0.0, 0.0, 1.0, 1.0]
+#domain.initialise_visualiser(rect=rect)
+domain.visualise = False
+
+
+domain.default_order = 2
+
+print domain.quantities.keys()
+
+domain.set_boundary( {'left': D, 'right': T, 'bottom': T, 'top': T} )
+domain.check_integrity()
+
+#Check that the boundary value gets propagated to all elements
+for t in domain.evolve(yieldstep = 0.1, finaltime = 1.5):
+    domain.write_time()
+    #pdb.set_trace()
Index: /trunk/anuga_core/demos/parallel/run_parallel_advection.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_advection.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_advection.py	(revision 9507)
@@ -0,0 +1,164 @@
+#!/usr/bin/env python
+#########################################################
+#
+#  Main file for parallel mesh testing. Runs an advection
+# flow simulation using a rectangular mesh
+#
+#
+#  Authors: Steve Roberts June 2005
+#  Modified by Linda Stals April 2006
+#
+#
+#########################################################
+
+
+
+# Parallel communication routines
+
+import pypar
+
+#  Mesh partition routines
+
+from parallel_meshes import parallel_rectangle
+
+# Parallel Domain
+ 
+from parallel_advection import Parallel_Domain
+from parallel_advection import Transmissive_boundary
+
+############################
+# Set the initial conditions
+############################
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+###############################
+# Read in processor information
+###############################
+
+numprocs = pypar.size()
+myid = pypar.rank()
+processor_name = pypar.get_processor_name()
+
+N = 5
+M = 2
+
+N = 10
+M = 5
+
+#######################
+# Partition the mesh
+#######################
+
+# Build a unit mesh, subdivide it over numproces processors with each
+# submesh containing M*N nodes
+
+points, vertices, boundary, full_send_dict, ghost_recv_dict =  \
+    parallel_rectangle(N, M, len1_g=1.0)
+
+#print 'ghost ',ghost_recv_dict
+#print 'full  ',full_send_dict
+
+
+
+print "Myid = ", myid, "no points = ", len(points), \
+      "no vertices = ", len(vertices), "no boundaries = ", len(boundary)
+
+###########################################
+# Start the computations on each subpartion
+###########################################
+
+# Create advection domain with direction (1,-1)
+# Initial condition is zero by default
+
+domain = Parallel_Domain(points, vertices, boundary,
+                         full_send_dict, ghost_recv_dict, velocity=[1.0, 0.0])
+
+#print 'ghost ',ghost_recv_dict
+#print 'full  ',full_send_dict
+
+
+# Boundaries
+
+T = Transmissive_boundary(domain)
+
+
+
+
+domain.set_boundary( {'left': T, 'right': T, 'bottom': T, 'top': T, \
+                      'ghost': None} )
+
+
+# Set Evolve parameters
+domain.set_default_order(2)
+domain.set_timestepping_method('rk2')
+
+print domain.get_timestepping_method()
+
+#domain.use_edge_limiter = True
+#domain.tight_slope_limiters = True
+#domain.use_centroid_velocities = False
+
+domain.CFL = 1.0
+
+domain.set_beta(0.8)
+
+
+# Ensure that the domain definitions make sense
+
+domain.check_integrity()
+
+# Set the inititial conditions
+
+domain.set_quantity('stage', Set_Stage(0.2,0.4,1.0))
+
+# Let processor 0 output some timing information
+
+visualise = False
+if visualise:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("stage", zScale = 5.0, dynamic=True)
+    vis.colour_height_quantity('stage', (0.2, 0.2, 0.8))
+    vis.start()
+    import time
+    time.sleep(2.0)
+    
+
+
+if myid == 0:
+    import time
+    t0 = time.time()
+
+for t in domain.evolve(yieldstep = 5.0, finaltime = 30.0):
+    if myid == 0:
+        domain.write_time()
+        
+    if visualise:
+        vis.update()						
+
+        
+
+if visualise: vis.evolveFinished()
+    
+
+# Output some computation statistics
+
+if myid == 0:
+    print 'That took %.2f seconds' %(time.time()-t0)
+    print 'Communication time %.2f seconds'\
+          %domain.communication_time
+    print 'Reduction Communication time %.2f seconds'\
+          %domain.communication_reduce_time
+        
+  
+if visualise: vis.join()
Index: /trunk/anuga_core/demos/parallel/run_parallel_advection_prof.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_advection_prof.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_advection_prof.py	(revision 9507)
@@ -0,0 +1,122 @@
+#!/usr/bin/env python
+#########################################################
+#
+#  Main file for parallel mesh testing. Runs an advection
+# flow simulation using a rectangular mesh
+#
+#
+#  Authors: Steve Roberts June 2005
+#  Modified by Linda Stals April 2006
+#
+#
+#########################################################
+
+
+
+# Parallel communication routines
+
+import pypar
+
+#  Mesh partition routines
+
+from parallel_meshes import parallel_rectangle
+
+# Parallel Domain
+ 
+from parallel_advection import Parallel_Domain
+from parallel_advection import Transmissive_boundary
+
+############################
+# Set the initial conditions
+############################
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+###############################
+# Read in processor information
+###############################
+
+numprocs = pypar.size()
+myid = pypar.rank()
+processor_name = pypar.Get_processor_name()
+
+N = 80
+M = 80
+
+#######################
+# Partition the mesh
+#######################
+
+# Build a unit mesh, subdivide it over numproces processors with each
+# submesh containing M*N nodes
+
+points, vertices, boundary, full_send_dict, ghost_recv_dict =  \
+    parallel_rectangle(N, M, len1_g=1.0)
+
+print "Myid = ", myid, "no points = ", len(points), \
+      "no vertices = ", len(vertices), "no boundaries = ", len(boundary)
+
+###########################################
+# Start the computations on each subpartion
+###########################################
+
+# Create advection domain with direction (1,-1)
+# Initial condition is zero by default
+
+domain = Parallel_Domain(points, vertices, boundary,
+                         full_send_dict, ghost_recv_dict, velocity=[1.0, 0.0])
+
+# Boundaries
+
+T = Transmissive_boundary(domain)
+domain.default_order = 2
+domain.set_boundary( {'left': T, 'right': T, 'bottom': T, 'top': T, \
+                      'ghost': None} )
+
+# Ensure that the domain definitions make sense
+
+domain.check_integrity()
+
+# Set the inititial conditions
+
+domain.set_quantity('stage', Set_Stage(0.2,0.4,1.0))
+
+# Let processor 0 output some timing information
+
+if myid == 0:
+    import time
+    t0 = time.time()
+
+# Start the evolve computions
+import hotshot
+profiler = hotshot.Profile("hotshot." + str(numprocs) + "." + str(myid) + ".prof")
+s = '''for t in domain.evolve(yieldstep = 0.1, finaltime = 2.0):
+  if myid == 0:
+    domain.write_time()
+'''
+    
+
+result = profiler.runctx(s, globals(), locals())
+profiler.close()
+
+# for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
+#     if myid == 0:
+#         domain.write_time()
+
+# Output some computation statistics
+
+if myid == 0:
+    print 'That took %.2f seconds' %(time.time()-t0)
+    print 'Communication time %.2f seconds'\
+          %domain.communication_time
+    print 'Reduction Communication time %.2f seconds'\
+          %domain.communication_reduce_time
Index: /trunk/anuga_core/demos/parallel/run_parallel_boyd_box_op.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_boyd_box_op.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_boyd_box_op.py	(revision 9507)
@@ -0,0 +1,167 @@
+import os.path
+import sys
+
+
+import anuga
+
+from anuga import myid, numprocs, finalize
+from anuga import Inlet_operator, Boyd_box_operator
+
+
+
+
+"""
+This test exercises the parallel culvert
+"""
+verbose = True
+
+    
+
+length = 40.
+width = 16.
+
+dx = dy = 2           # Resolution: Length of subdivisions on both axes
+
+#----------------------------------------------------------------------
+# Setup initial conditions
+#----------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=-x/1000
+    
+    N = len(x)
+    for i in range(N):
+
+       # Sloping Embankment Across Channel
+        if 5.0 < x[i] < 10.1:
+            # Cut Out Segment for Culvert face                
+            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+               z[i]=z[i]
+            else:
+               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+        if 10.0 < x[i] < 12.1:
+           z[i] +=  2.5                    # Flat Crest of Embankment
+        if 12.0 < x[i] < 14.5:
+            # Cut Out Segment for Culvert face                
+            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+               z[i]=z[i]
+            else:
+               z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                   
+        
+    return z
+
+##-----------------------------------------------------------------------
+## Setup domain
+##-----------------------------------------------------------------------
+if myid == 0:
+    points, vertices, boundary = anuga.rectangular_cross(int(length/dx),
+                                                   int(width/dy),
+                                                   len1=length, 
+                                                   len2=width)
+
+    domain = anuga.Domain(points, vertices, boundary)   
+    domain.set_name()                 # Output name output_script_name.sww
+    domain.set_flow_algorithm('1_5')
+else:
+    domain = None
+    
+##-----------------------------------------------------------------------
+## Distribute domain
+##-----------------------------------------------------------------------
+
+
+domain = anuga.distribute(domain)
+#domain.dump_triangulation("run_parallel_boyd_box_op_domain.png")
+    
+
+##-----------------------------------------------------------------------
+## Setup boundary conditions
+##-----------------------------------------------------------------------
+    
+domain.set_quantity('elevation', topography) 
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                        expression='elevation')   # Dry initial condition
+
+    
+Bi = anuga.Dirichlet_boundary([5.0, 0.0, 0.0])
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+
+    ################ Define Fractional Operators ##########################
+line0 = [[10.0, 10.0], [30.0, 10.0]]
+#line0 = [[29.0, 10.0], [30.0, 10.0]]
+poly1 = [[0.0, 10.0], [0.0, 15.0], [5.0, 15.0], [5.0, 10.0]]
+Q0 = anuga.file_function('../data/test_hydrograph.tms', quantities=['hydrograph'])
+Q1 = 5.0
+
+samples = 50
+
+
+inlet0 = None
+inlet1 = None
+boyd_box0 = None
+
+inlet0 = Inlet_operator(domain, line0, Q0, logging=True, description='inlet0', verbose = False)
+inlet1 = Inlet_operator(domain, poly1, Q1, logging=True, description='inlet1', verbose = False)
+
+# Enquiry point [ 19.    2.5] is contained in two domains in 4 proc case
+
+boyd_box0 = Boyd_box_operator(domain,
+                              end_points=[[9.0, 2.5],[19.0, 2.5]],
+                              losses=1.5,
+                              width=5.0,
+                              apron=5.0,
+                              use_momentum_jet=True,
+                              use_velocity_head=False,
+                              manning=0.013,
+                              logging=True,
+                              description='boyd_box_0',
+                              verbose=False)
+
+if inlet0 is not None and verbose: inlet0.print_statistics()
+if inlet1 is not None and verbose: inlet1.print_statistics()
+if boyd_box0 is not None and verbose: boyd_box0.print_statistics()
+
+
+##-----------------------------------------------------------------------
+## Evolve system through time
+##-----------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 2.0, finaltime = 20.0):
+    if myid == 0 and verbose:
+        domain.write_time()
+
+    #print domain.volumetric_balance_statistics()
+
+    stage = domain.get_quantity('stage')
+
+
+    if boyd_box0 is not None and verbose : boyd_box0.print_timestepping_statistics()
+
+    #for i in range(samples):
+    #    if tri_ids[i] >= 0:                
+    #        if verbose: print 'P%d tri %d, value = %s' %(myid, i, stage.centroid_values[tri_ids[i]])
+
+    sys.stdout.flush()
+
+    pass
+
+##-----------------------------------------------------------------------
+## Assign/Test Control data
+##-----------------------------------------------------------------------
+
+domain.sww_merge(delete_old=True)
+
+finalize()
+
+    
+
Index: /trunk/anuga_core/demos/parallel/run_parallel_gate_operator.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_gate_operator.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_gate_operator.py	(revision 9507)
@@ -0,0 +1,192 @@
+import os.path
+import sys
+
+#from anuga.utilities.system_tools import get_pathname_from_package
+#from anuga.geometry.polygon_function import Polygon_function
+        
+#from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+#from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+import anuga
+import time
+
+                            
+from math import pi, pow, sqrt
+
+import numpy as num
+
+from anuga import distribute, myid, numprocs, finalize, barrier
+
+from anuga import Inlet_operator, Boyd_box_operator
+
+
+
+#----------------------------------------------------------------------
+# Setup initial conditions
+#----------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=-x/1000
+
+    N = len(x)
+    for i in range(N):
+
+        # Sloping Embankment Across Channel
+        if 5.0 < x[i] < 10.1:
+            # Cut Out Segment for Culvert face
+            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0:
+                z[i]=z[i]
+            else:
+                z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+        if 10.0 < x[i] < 12.1:
+            z[i] +=  2.5                    # Flat Crest of Embankment
+        if 12.0 < x[i] < 14.5:
+            # Cut Out Segment for Culvert face
+            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+                z[i]=z[i]
+            else:
+                z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+
+    return z
+
+  
+
+length = 40.
+width = 15.
+
+dx = dy = 0.5          # Resolution: Length of subdivisions on both axes
+
+
+if myid == 0:
+
+    points, vertices, boundary = anuga.rectangular_cross(int(length/dx),
+                                                         int(width/dy),
+                                                         len1=length, 
+                                                         len2=width)
+    domain = anuga.Domain(points, vertices, boundary)
+    domain.set_name()                 # Output name
+    domain.set_flow_algorithm('DE0')
+    domain.set_store_vertices_uniquely(True)
+    domain.set_quantity('elevation', topography)
+    domain.set_quantity('friction', 0.01)         # Constant friction
+    domain.set_quantity('stage',
+                        expression='elevation')   # Dry initial condition
+
+
+
+else:
+
+    domain = None
+
+
+domain = distribute(domain, verbose = True)
+
+#domain.set_store_vertices_uniquely(False)
+
+
+gate = Boyd_box_operator(domain,
+                            end_points=[[9.0, 2.5],[13.0, 2.5]],
+                            losses=1.5,
+                            width=1.5,
+                            height = 0.0001,
+                            apron=5.0,
+                            use_momentum_jet=True,
+                            use_velocity_head=False,
+                            manning=0.013,
+                            verbose=False)
+
+# Close the gate
+if gate is not None:
+    gate.set_culvert_height(0.0)
+
+line = [[0.0, 5.0], [0.0, 10.0]]
+Q = 1.0
+in0 = Inlet_operator(domain, line, Q)
+
+
+
+
+##-----------------------------------------------------------------------
+## Setup boundary conditions
+##-----------------------------------------------------------------------
+
+## Inflow based on Flow Depth and Approaching Momentum
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+##-----------------------------------------------------------------------
+## Evolve system through time
+##-----------------------------------------------------------------------
+barrier()
+
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 1.0, finaltime = 50):
+
+    if myid == 0:
+        print 80*'='
+        domain.write_time()
+
+    #================================================
+    # gate operator lives on a number of processors
+    #
+    # Need to run parallel gate operations on all
+    # processors where gate lives.
+    # But true result only on master processor.
+    #================================================
+    if gate is not None:
+        [s0, s1] = gate.get_enquiry_stages()
+        [d0, d1] = gate.get_enquiry_depths()
+        [e0, e1] = gate.get_enquiry_elevations()
+        [i0, i1] = gate.get_enquiry_invert_elevations()
+        [w0, w1] = gate.get_enquiry_water_depths()
+
+        output = gate.discharge_routine()
+
+        if myid == gate.get_master_proc():
+            print 'myid ', myid, s0,s1
+            print 'myid ', myid, d0,d1
+            print 'myid ', myid, e0,e1
+            print 'myid ', myid, i0,i1
+            print 'myid ', myid, w0,w1
+
+            print 'myid ',myid, output
+
+
+            if d0 > 0.2: gate.set_culvert_height(10.0)
+
+            
+    
+
+barrier()
+
+
+domain.sww_merge(delete_old=True)
+
+
+
+
+barrier()
+
+for p in range(numprocs):
+    if myid == p:
+        print 'Processor %g ' %myid
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'Communication time %.2f seconds'%domain.communication_time
+        print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+        print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+    else:
+        pass
+
+    barrier()
+
+
+finalize()
+
Index: /trunk/anuga_core/demos/parallel/run_parallel_job.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_job.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_job.py	(revision 9507)
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+
+"""Run job in parallel using PBS batch system with
+   specified number of processors
+
+   Syntax
+     run_parallel_job.py 4 <prog>
+"""
+
+queue = 'normal'
+walltime = '00:10:00'
+vmem = '1Mb'
+
+import sys, string, time, os
+arg = sys.argv
+command = arg[0]
+ext = '.py'
+
+if len(arg)<3:
+    print __doc__
+    print 'Usage:'
+    print command + ' p <filename>' + ext
+    sys.exit()
+
+numproc = int(arg[1])    
+progname = arg[2]
+
+
+#Create batchfile
+
+s = """
+#!/bin/csh
+
+# This file is automatically generated - don't edit
+#
+
+#PBS -wd  
+#PBS -q %s
+#PBS -l ncpus=%d
+#PBS -l walltime=%s,vmem=%s
+#PBS -l software=python
+
+prun %s""" %(queue, numproc, walltime, vmem, progname)
+
+
+scriptname = progname + '.run'
+fid = open(scriptname, 'w')
+fid.write(s)
+fid.close()
+
+#Cleanup
+os.system('/bin/rm *.o* *.e*')
+
+# Run
+os.system('qsub %s' %scriptname)
+os.system('/bin/rm %s' %scriptname)
+
+#Create check scripts
+fid = open('q', 'w')
+fid.write("""#!/bin/csh
+#This file is generated automatically.
+#Any modifications will be lost.
+#Edit runjob.py instead.
+echo '--------------------------------------------'
+nqstat %s\n""" %queue)
+fid.close()
+
+#fid = open('o', 'w')
+#fid.write("""#!/bin/csh
+##This file is generated automatically.
+##Any modifications will be lost.
+##Edit runjob.py instead.
+#cat %s.e* %s.o*\n""" %((progname,)*2))
+#fid.close()
+
+
+fid = open('o', 'w')
+fid.write("""#!/bin/csh
+#This file is generated automatically.
+#Any modifications will be lost.
+#Edit runjob.py instead.
+cat *.e* *.o*\n""")
+
+fid.close()
+
+os.system('chmod +x o q')
+
+
+
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/parallel/run_parallel_merimbula_test.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_merimbula_test.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_merimbula_test.py	(revision 9507)
@@ -0,0 +1,173 @@
+#!/usr/bin/env python
+#########################################################
+#
+#  Main file for parallel mesh testing. Runs an advection
+# flow simulation using a rectangular mesh
+#
+#  This is a modification of the run_parallel_advection.py
+# file
+#
+#
+# *) The (new) files that have been added to manage the
+# grid partitioning are
+#    +) pmesh_divide_metis.py: subdivide a pmesh
+#    +) build_submesh.py: build the submeshes on the host
+# processor.
+#    +) build_local.py: build the GA mesh datastructure
+# on each processor.
+#    +) build_commun.py: handle the communication between
+# the host and processors
+#
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+
+import pypar    # The Python-MPI interface
+import time
+
+from sys import path
+
+print path
+
+# Numeric arrays
+from Numeric import array, zeros, Float
+
+from print_stats import print_test_stats, build_full_flag
+
+from anuga.abstract_2d_finite_volumes.pmesh2domain\
+     import pmesh_to_domain_instance
+from anuga.advection.advection import Domain as Advection_Domain
+from parallel_advection import Parallel_Domain
+
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Transmissive_boundary
+
+# mesh partition routines
+
+from pmesh_divide  import pmesh_divide_metis
+from build_submesh import build_submesh
+from build_local   import build_local_mesh
+from build_commun  import send_submesh, rec_submesh, extract_submesh
+
+
+
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+# read in the processor information
+
+numprocs = pypar.size()
+myid = pypar.rank()
+processor_name = pypar.Get_processor_name()
+
+
+#-------
+# Domain
+rect = zeros( 4, Float) # Buffer for results
+
+if myid == 0:
+
+    # read in the test files
+
+    filename = 'test-100.tsh'
+#    filename = 'merimbula_10785.tsh'
+    nx = numprocs
+    ny = 1
+    if nx*ny != numprocs:
+        print "WARNING: number of subboxes is not equal to the number of proc"
+
+    domain_full = pmesh_to_domain_instance(filename, Advection_Domain)
+    domain_full.set_quantity('stage', Set_Stage(200.0,300.0,1.0))
+#    domain_full.set_quantity('stage', Set_Stage(756000.0,756500.0,4.0))
+
+    nodes, triangles, boundary, triangles_per_proc, quantities  =\
+            pmesh_divide_metis(domain_full, numprocs)
+    
+    # subdivide the mesh
+    rect = array(domain_full.xy_extent, Float)
+#    rect = array(rect, Float)
+
+    submesh = build_submesh(nodes, triangles, boundary, quantities, \
+                            triangles_per_proc)
+
+    # send the mesh partition to the appropriate processor
+
+    for p in range(1, numprocs):
+      send_submesh(submesh, triangles_per_proc, p)
+
+    # Build the local mesh for processor 0
+    
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
+              extract_submesh(submesh, triangles_per_proc)
+
+# read in the mesh partition that belongs to this
+# processor (note that the information is in the
+# correct form for the GA data structure
+
+else:
+    [points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict] = \
+             rec_submesh(0)
+
+pypar.broadcast(rect,0)
+
+domain = Parallel_Domain(points, vertices, boundary,
+                                   full_send_dict  = full_send_dict,
+                                   ghost_recv_dict = ghost_recv_dict,
+                                   velocity = [0.1,0.0])
+
+# Make a notes of which triangles are full and which are ghost
+
+tri_full_flag = build_full_flag(domain, ghost_recv_dict)
+
+#domain.initialise_visualiser(rect=rect)
+
+#Boundaries
+
+T = Transmissive_boundary(domain)
+#R = Reflective_boundary(domain)
+domain.set_boundary( {'outflow': T, 'inflow': T, 'inner':T, 'exterior': T, 'open':T, 'ghost':None} )
+
+
+
+domain.set_quantity('stage', quantities['stage'])
+
+#---------
+# Evolution
+t0 = time.time()
+
+
+try:
+    domain.initialise_visualiser(rect=rect)
+    #domain.visualiser.coloring['stage'] = True
+    #domain.visualiser.scale_z['stage'] = 0.2
+except:
+    print 'No visualiser'
+
+
+from norms import linf_norm
+
+yieldstep = 1
+finaltime = 200
+
+#yieldstep = 1000
+#finaltime = 50000
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+    print_test_stats(domain, tri_full_flag)
+
+if myid == 0:
+    print 'That took %.2f seconds' %(time.time()-t0)
Index: /trunk/anuga_core/demos/parallel/run_parallel_sw_merimbula.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_sw_merimbula.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_sw_merimbula.py	(revision 9507)
@@ -0,0 +1,200 @@
+
+"""Run parallel shallow water domain.
+
+   run using command like:
+
+   mpirun -np m python run_parallel_sw_merimbula.py
+
+   where m is the number of processors to be used.
+   
+   Will produce sww files with names domain_Pn_m.sww where m is number of processors and
+   n in [0, m-1] refers to specific processor that owned this part of the partitioned mesh.
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import os
+import sys
+import time
+import numpy as num
+
+#------------------------
+# ANUGA Modules
+#------------------------
+	
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+from anuga import Transmissive_boundary
+from anuga import Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
+
+from anuga import rectangular_cross
+from anuga import create_domain_from_file
+
+
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
+
+
+#--------------------------------------------------------------------------
+# Setup parameters
+#--------------------------------------------------------------------------
+
+mesh_filename = "data/merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0; yieldstep = 10; finaltime = 100
+#mesh_filename = "data/merimbula_17156.tsh"   ; x0 = 756000.0 ; x1 = 756500.0; yieldstep = 50; finaltime = 500
+#mesh_filename = "data/merimbula_43200_1.tsh"   ; x0 = 756000.0 ; x1 = 756500.0; yieldstep = 50; finaltime = 500
+#mesh_filename = "data/test-100.tsh" ; x0 = 200.0 ; x1 = 300.0; yieldstep = 1; finaltime = 10
+#mesh_filename = "data/test-20.tsh" ; x0 = 250.0 ; x1 = 350.0; yieldstep = 1; finaltime = 50
+
+verbose = False
+
+#--------------------------------------------------------------------------
+# Setup procedures
+#--------------------------------------------------------------------------
+class Set_Stage:
+    """Set an initial condition with constant water height, for x0<x<x1
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))+1.0
+
+
+class Set_Elevation:
+    """Set an elevation
+    """
+
+    def __init__(self, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return x/self.h
+    
+
+#--------------------------------------------------------------------------
+# Setup Domain only on processor 0
+#--------------------------------------------------------------------------
+if myid == 0:
+    domain = create_domain_from_file(mesh_filename)
+    domain.set_quantity('stage', Set_Stage(x0, x1, 1.0))
+    #domain.set_quantity('stage', 0.0)
+    #domain.set_datadir('.')
+    domain.set_name()
+    domain.set_store(True)
+    #domain.set_store_vertices_smoothly(True)
+    domain.set_store_vertices_smoothly(False)
+    #domain.set_quantity('elevation', Set_Elevation(500.0))
+
+    #print domain.statistics()
+    #print domain.get_extent()
+    #print domain.get_extent(absolute=True)
+    #print domain.geo_reference
+else:
+    domain = None
+
+#--------------------------------------------------------------------------
+# Distribute sequential domain on processor 0 to other processors
+#--------------------------------------------------------------------------
+
+if myid == 0 and verbose: print 'DISTRIBUTING DOMAIN'
+domain = distribute(domain, verbose=verbose)
+
+#--------------------------------------------------------------------------
+# On all processors, setup evolve parameters for domains on all processors
+# (all called "domain"
+#--------------------------------------------------------------------------
+
+domain.set_flow_algorithm('2_0')
+domain.set_store_centroids()
+
+#for p in range(numprocs):
+#    if myid == p:
+#        print 'P%d'%p
+#        print domain.get_extent()
+#        print domain.get_extent(absolute=True)
+#        print domain.geo_reference
+#        print domain.s2p_map
+#        print domain.p2s_map
+#        print domain.tri_l2g
+#        print domain.node_l2g
+#    else:
+#        pass
+#
+#    barrier()
+
+
+domain.set_quantities_to_be_stored({'elevation':1,
+                                    'friction':1,
+                                    'stage':2,
+                                    'xmomentum':2,
+                                    'ymomentum':2})
+                                 
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+# This must currently happen *after* domain has been distributed
+#------------------------------------------------------------------------------
+Br = Reflective_boundary(domain)      # Solid reflective wall
+
+from math import sin
+Bts = Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, lambda t: 10*sin(t/60))
+
+domain.set_boundary({'outflow' :Br, 'inflow' :Br, 'inner' :Br, 'exterior' :Br, 'open' :Bts})
+
+
+#bdy_ids = domain.tag_boundary_cells['open']
+#vol_ids = domain.boundary_cells[bdy_ids]
+#edge_ids = domain.boundary_edges[bdy_ids]
+
+#print domain.mesh.edge_midpoint_coordinates[3*vol_ids+edge_ids]
+
+
+
+#------------------------------------------------------------------------------
+# Evolution
+#------------------------------------------------------------------------------
+if myid == 0 and verbose: print 'EVOLVE'
+
+barrier()
+t0 = time.time()
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    #print 'P:'+ str(myid) + ' '+ domain.timestepping_statistics()
+    if myid == 0:
+        domain.write_time()
+
+
+barrier()
+
+for p in range(numprocs):
+    if myid == p:
+        print 'Processor %g ' %myid
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'Communication time %.2f seconds'%domain.communication_time
+        print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+        print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+    else:
+        pass
+
+    barrier()
+
+
+#--------------------------------------------------
+# Merge the individual sww files into one file
+#--------------------------------------------------
+domain.sww_merge(delete_old=True)
+
+
+domain.dump_triangulation()
+
+finalize()
+
+
+
+
Index: /trunk/anuga_core/demos/parallel/run_parallel_sw_merimbula_profile.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_sw_merimbula_profile.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_sw_merimbula_profile.py	(revision 9507)
@@ -0,0 +1,186 @@
+
+"""Run parallel shallow water domain.
+
+   run using command like:
+
+   mpirun -np m python run_parallel_sw_merimbula.py
+
+   where m is the number of processors to be used.
+   
+   Will produce sww files with names domain_Pn_m.sww where n is number of processors and
+   m in [0, n-1] refers to specific processor that owned this part of the partitioned mesh.
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import os
+import sys
+import time
+import numpy as num
+
+#------------------------
+# ANUGA Modules
+#------------------------
+	
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+from anuga import Transmissive_boundary
+
+from anuga import rectangular_cross
+from anuga import create_domain_from_file
+
+
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
+
+
+#--------------------------------------------------------------------------
+# Setup parameters
+#--------------------------------------------------------------------------
+
+#mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0
+mesh_filename = "merimbula_43200.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
+#mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
+
+finaltime = 500
+yieldstep = 50
+verbose = True
+
+#--------------------------------------------------------------------------
+# Setup procedures
+#--------------------------------------------------------------------------
+class Set_Stage:
+    """Set an initial condition with constant water height, for x0<x<x1
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1)) + 1.0
+
+#--------------------------------------------------------------------------
+# Setup Domain only on processor 0
+#--------------------------------------------------------------------------
+if myid == 0:
+    domain = create_domain_from_file(mesh_filename)
+    domain.set_quantity('stage', Set_Stage(x0, x1, 2.0))
+else:
+    domain = None
+
+#--------------------------------------------------------------------------
+# Distribute sequential domain on processor 0 to other processors
+#--------------------------------------------------------------------------
+
+if myid == 0 and verbose: print 'DISTRIBUTING DOMAIN'
+domain = distribute(domain)
+
+#domain.smooth = False
+barrier()
+for p in range(numprocs):
+    if myid == p:
+        print 'Process ID %g' %myid
+        print 'Number of triangles %g ' %domain.get_number_of_triangles()
+
+    barrier()
+
+
+domain.set_flow_algorithm(2.0)
+
+if myid == 0:
+    domain.print_algorithm_parameters()
+
+barrier()
+
+domain.set_name('meribula')
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+# This must currently happen *after* domain has been distributed
+#------------------------------------------------------------------------------
+Br = Reflective_boundary(domain)      # Solid reflective wall
+
+domain.set_boundary({'exterior' :Br, 'open' :Br})
+
+
+barrier()
+for p in range(numprocs):
+    if myid == p:
+        print domain.boundary_statistics()
+
+    barrier()
+
+#------------------------------------------------------------------------------
+# Evolution
+#------------------------------------------------------------------------------
+if myid == 0 and verbose: print 'EVOLVE'
+
+t0 = time.time()
+
+
+s = """
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+"""
+
+
+# Profiling
+import cProfile
+prof_file = 'evolve-prof'+ str(numprocs) + '_' + str(myid) + '.dat'
+cProfile.run(s,prof_file)
+
+
+barrier()
+
+if myid == 0:
+    import pstats
+    p = pstats.Stats(prof_file)
+    #p.sort_stats('cumulative').print_stats(25)
+
+
+    p.sort_stats('time').print_stats(25)
+
+barrier()
+
+
+if myid == 1:
+    import pstats
+    p = pstats.Stats(prof_file)
+    #p.sort_stats('cumulative').print_stats(25)
+
+
+    p.sort_stats('time').print_stats(25)
+
+    #p.print_stats()
+
+# Run evolve loop
+#result = profiler.runctx(s, globals(), locals())
+
+#result.dump_stats("profile." + str(numprocs) + "." + str(myid) + ".dat")
+#profiler.close()
+
+barrier()
+
+for p in range(numprocs):
+    if myid == p:
+        print 'Process ID %g' %myid
+        print 'Number of processors %g ' %numprocs
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'Communication time %.2f seconds'%domain.communication_time
+        print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+        print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+
+    barrier()
+
+
+finalize()
+
+
+
+
Index: /trunk/anuga_core/demos/parallel/run_parallel_sw_merimbula_test.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_sw_merimbula_test.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_sw_merimbula_test.py	(revision 9507)
@@ -0,0 +1,283 @@
+#!/usr/bin/env python
+#########################################################
+#
+#  Main file for parallel mesh testing.
+#
+#  This is a modification of the run_parallel_advection.py
+# file.
+#
+#
+# *) The (new) files that have been added to manage the
+# grid partitioning are
+#    +) pmesh_divide_metis.py: subdivide a pmesh
+#    +) build_submesh.py: build the submeshes on the host
+# processor.
+#    +) build_local.py: build the GA mesh datastructure
+# on each processor.
+#    +) build_commun.py: handle the communication between
+# the host and processors
+#
+# *) Things still to do:
+#    +) Overlap the communication and computation: The
+# communication routines in build_commun.py should be
+# interdispersed in the build_submesh.py and build_local.py
+# files. This will overlap the communication and
+# computation and will be far more efficient. This should
+# be done after more testing and there more confidence in
+# the subpartioning.
+#    +) Much more testing especially with large numbers of
+# processors.
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+import sys
+import pypar    # The Python-MPI interface
+import time
+
+# Numeric arrays
+import numpy as num
+#from numpy import array, zeros, float
+
+# Print debugging information
+from print_stats import print_test_stats, build_full_flag
+
+# pmesh
+from anuga.shallow_water import Domain
+from parallel_shallow_water import Parallel_domain
+from anuga.abstract_2d_finite_volumes.pmesh2domain\
+     import pmesh_to_domain_instance
+
+# Reuse previous mesh import
+from anuga.caching import cache
+
+# Mesh partition routines
+from distribute_mesh  import pmesh_divide_metis
+from distribute_mesh  import build_submesh
+from distribute_mesh  import build_local_mesh
+from distribute_mesh  import send_submesh, rec_submesh, extract_submesh
+
+
+###############################
+# Read in processor information
+###############################
+
+numprocs = pypar.size()
+myid = pypar.rank()
+processor_name = pypar.get_processor_name()
+
+############################
+# Set the initial conditions
+############################
+
+rect = num.zeros( 4, num.float) # Buffer for results
+
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+#######################
+# Partition the domain
+#######################
+
+if myid == 0:
+
+    # Read in the test files
+
+#    filename = 'test-100.tsh'
+#    filename = 'merimbula_10785_1.tsh'
+    filename = 'merimbula_43200.tsh'
+
+    # Build the whole domain
+    
+    domain_full = pmesh_to_domain_instance(filename, Domain)
+
+#    domain_full = cache(pmesh_to_domain_instance,
+#               (filename, Domain),
+#              dependencies = [filename])
+
+    rect = num.array(domain_full.get_extent(), num.float)
+    print rect
+
+    # Initialise the wave
+
+    #domain_full.set_quantity('stage', Set_Stage(200.0,300.0,1.0))
+    domain_full.set_quantity('stage', Set_Stage(756000.0,756500.0,2.0))
+#    domain_full.set_quantity('stage', Set_Stage(756000.0,756500.0,0.0))
+
+    # Subdivide the domain
+
+    # Note the different arguments compared with pmesh_divide,
+    # pmesh_divide_steve etc.
+    
+    nodes, triangles, boundary, triangles_per_proc, quantities = \
+         pmesh_divide_metis(domain_full, numprocs)
+
+    print triangles_per_proc
+    
+    rect = num.array(domain_full.get_extent(), num.float)
+
+    submesh = build_submesh(nodes, triangles, boundary,\
+                            quantities, triangles_per_proc)
+
+    # Send the mesh partition to the appropriate processor
+
+    for p in range(1, numprocs):
+      send_submesh(submesh, triangles_per_proc, p)
+
+    # Build the local mesh for processor 0
+
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
+             extract_submesh(submesh, triangles_per_proc)
+
+# Read in the mesh partition that belongs to this
+# processor (note that the information is in the
+# correct form for the GA data structure
+
+else:
+    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict , \
+            no_full_nodes, no_full_trigs = rec_submesh(0)
+
+
+###########################################
+# Start the computations on each subpartion
+###########################################
+
+#if myid == 0:
+#    print 'ghost'
+#    print ghost_recv_dict
+#processor_name
+#if myid == 0:
+#    print 'full'
+#    print full_send_dict
+
+# The visualiser needs to know the size of the whole domain
+
+pypar.broadcast(rect,0)
+
+domain = Parallel_domain(points, vertices, boundary,
+                         full_send_dict  = full_send_dict,
+                         ghost_recv_dict = ghost_recv_dict)
+
+# Make a note of which triangles are full and which are ghost
+
+tri_full_flag = build_full_flag(domain, ghost_recv_dict)
+
+try:
+    #domain.initialise_visualiser(rect=rect)
+    #domain.visualiser.coloring['stage'] = True
+    #domain.visualiser.scale_z['stage'] = 0.2
+    #domain.visualiser.scale_z['elevation'] = 0.05
+    pass
+except:
+    print 'No visualiser'
+
+
+
+domain.default_order = 1
+
+#Boundaries
+from anuga.interface import Transmissive_boundary, Reflective_boundary
+
+T = Transmissive_boundary(domain)
+R = Reflective_boundary(domain)
+domain.set_boundary( {'outflow': R, 'inflow': R, 'inner':R, 'exterior': R, 'open':R, 'ghost':None} )
+
+
+domain.set_quantity('stage', quantities['stage'])
+domain.set_quantity('elevation', quantities['elevation'])
+
+domain.store = False
+
+#---------
+# Evolution
+t0 = time.time()
+
+print 'Processor %d on %s: No of elements %d'%(domain.processor,processor_name,domain.number_of_elements)
+yieldstep = 50.0
+finaltime = 500.0
+
+#yieldstep = 1000
+#finaltime = 40000
+
+#yieldstep = 1
+#finaltime = 1
+#processor_name
+#for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+#    if myid == 0:
+#        domain.write_time()
+        #print 'Processor %d, Integral of stage %d'%\
+        #       (domain.processor,domain.quantities['stage'].get_integral())
+        #    print_test_stats(domain, tri_full_flag)
+
+
+# Profiling
+#import profile
+#profiler = profile.Profile()
+#result.dump_stats("profile." + str(numprocs) + "." + str(myid) + ".dat")
+
+## #New hotshot profiling
+## import hotshot
+## profiler = hotshot.Profile("hotshot." + str(numprocs) + "." + str(myid) + ".prof")
+## s = '''for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+##   if myid == 0:
+##     domain.write_time()
+##   print_test_stats(domain, tri_full_flag)
+
+## '''
+## result = profiler.runctx(s, globals(), locals())
+## profiler.close()
+
+#from vtk_realtime_visualiser import Visualiser
+#V = Visualiser(domain,default_scale_z=100.0)
+#V.coloring['stage'] = True
+#V.coloring['elevation'] = False
+#V.setup['elevation']=True
+#V.updating['stage']=True
+#V.qcolor['stage'] = (0.1,0.4,0.99)
+
+
+#V.start()
+#V.idle.wait()
+#V.idle.clear()
+
+
+
+
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+        #print_test_stats(domain, tri_full_flag)
+
+#    V.redraw_ready.set()
+#    V.idle.wait()
+#    V.idle.clear()
+#    V.unpaused.wait()
+
+
+#print 'P%d: That took %.2f seconds' %(myid, time.time()-t0)
+#print 'P%d: Communication time %.2f seconds' %(myid, domain.communication_time)
+#print 'P%d: Reduction Communication time %.2f seconds' %(myid, domain.communication_reduce_time)
+#print 'P%d: Broadcast time %.2f seconds' %(myid, domain.communication_broadcast_time)
+
+
+
+if myid == 0:
+    print 'That took %.2f seconds' %(time.time()-t0)
+    print 'Communication time %.2f seconds'%domain.communication_time
+    print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+    print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+
+
+pypar.finalize()
Index: /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangle.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangle.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangle.py	(revision 9507)
@@ -0,0 +1,162 @@
+#########################################################
+#
+#  Main file for parallel mesh testing.
+#
+#  This is a modification of the run_parallel_advection.py
+# file.
+#
+#
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+
+import time
+
+import numpy as num
+
+from print_stats import print_test_stats, build_full_flag
+
+#----------------------------
+# Sequential interface
+#---------------------------
+from anuga import Domain
+from anuga import Transmissive_boundary, Reflective_boundary
+
+#----------------------------
+# Parallel interface
+#---------------------------
+from anuga_parallel import Parallel_shallow_water_domain
+from anuga_parallel import parallel_rectangle
+from anuga_parallel import myid, numprocs, finalize, get_processor_name
+
+
+processor_name = get_processor_name()
+
+M = 50
+N = M*numprocs
+
+
+if myid == 0:
+    print 'N == %d' %N
+
+points, vertices, boundary, full_send_dict, ghost_recv_dict =\
+        parallel_rectangle(N, M, len1_g=1.0*numprocs, len2_g = 1.0)
+
+print "Myid = ", myid, "no points = ", len(points), \
+      "no vertices = ", len(vertices), "no boundaries = ", len(boundary)
+
+###########################################
+# Start the computations on each subpartion
+###########################################
+
+domain = Parallel_shallow_water_domain(points, vertices, boundary,
+                                       full_send_dict  = full_send_dict,
+                                       ghost_recv_dict = ghost_recv_dict)
+
+
+print 'after parallel domain'
+
+
+
+domain.set_name('sw_rectangle')
+
+#Boundaries
+T = Transmissive_boundary(domain)
+R = Reflective_boundary(domain)
+
+
+domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R, 'ghost': None} )
+
+
+print 'after set_boundary'
+
+
+
+domain.check_integrity()
+
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.75, y0=0.0, y1=1.0, h=5.0, h0=0.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.y0 = y0
+        self.y1 = y1
+        self.h  = h
+        self.h0 = h0
+
+    def __call__(self, x, y):
+        return self.h0 + self.h*((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1))
+
+print 'after check_integrity'
+
+domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
+
+
+print 'after set quantity'
+
+# Set Evolve parameters
+domain.set_default_order(2)
+domain.set_timestepping_method('rk2')
+domain.set_CFL(0.7)
+domain.set_beta(1.5)
+
+#print domain.get_timestepping_method()
+#domain.use_edge_limiter = True
+#domain.tight_slope_limiters = True
+#domain.use_centroid_velocities = False
+
+
+print 'after evolve parameters'
+
+
+import pdb; pdb.set_trace()
+
+if myid == 0:
+    import time
+    t0 = time.time()
+
+
+# Turn on the visualisation
+visualise = False
+if visualise:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("elevation", offset=0.001, dynamic=False)
+    vis.render_quantity_height("stage", dynamic=True)
+    vis.colour_height_quantity('stage', (0.2, 0.2, 0.8))
+    vis.start()
+    import time
+    time.sleep(2.0)
+
+
+
+yieldstep = 0.05
+finaltime = 2.0
+
+#Check that the boundary value gets propagated to all elements
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+    #print_test_stats(domain, tri_full_flag)
+    if visualise:
+        vis.update()						
+
+
+
+if visualise: vis.evolveFinished()
+
+if myid == 0:
+    print 'That took %.2f seconds' %(time.time()-t0)
+    print 'Communication time %.2f seconds'%domain.communication_time
+    print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+    print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+
+
+if visualise: vis.join()
+
+finalize()
Index: /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangular_cross.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangular_cross.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangular_cross.py	(revision 9507)
@@ -0,0 +1,162 @@
+#########################################################
+#
+#  Main file for parallel mesh testing.
+#
+#  This is a modification of the run_parallel_advection.py
+# file.
+#
+#
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+
+import time
+import sys
+
+
+#----------------------------
+# Sequential interface
+#---------------------------
+from anuga import Transmissive_boundary, Reflective_boundary
+from anuga import rectangular_cross_domain
+from anuga import Set_quantity
+
+#----------------------------
+# Parallel interface
+#---------------------------
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
+
+
+t0 = time.time()
+
+verbose = True
+
+#--------------------------------------------------------------------------
+# Setup Domain only on processor 0
+#--------------------------------------------------------------------------
+if myid == 0:
+    length = 2.0
+    width = 2.0
+    #dx = dy = 0.005
+    #dx = dy = 0.00125
+    dx = dy  = 0.5
+    domain = rectangular_cross_domain(int(length/dx), int(width/dy),
+                                              len1=length, len2=width, origin=(-length/2, -width/2), verbose=verbose)
+
+
+    domain.set_store(True)
+    domain.set_quantity('elevation', lambda x,y : -1.0-x )
+    domain.set_quantity('stage', 1.0)
+    domain.set_flow_algorithm('DE0')
+    domain.print_statistics()
+else:
+    domain = None
+
+t1 = time.time()
+
+if myid == 0 :
+    print 'Create sequential domain: Time',t1-t0
+
+if myid == 0 and verbose: 
+    print 'DISTRIBUTING DOMAIN'
+    sys.stdout.flush()
+    
+barrier()
+
+#-------------------------------------------------------------------------
+# Distribute domain
+#-------------------------------------------------------------------------
+domain = distribute(domain,verbose=verbose)
+
+
+t2 = time.time()
+
+if myid == 0 :
+    print 'Distribute domain: Time ',t2-t1
+    
+if myid == 0 : print 'after parallel domain'
+
+
+
+domain.set_name('sw_rectangle')
+
+#Boundaries
+T = Transmissive_boundary(domain)
+R = Reflective_boundary(domain)
+
+
+domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R, 'ghost': None} )
+
+
+if myid == 0 : print 'after set_boundary'
+
+
+
+#domain.check_integrity()
+
+if myid == 0 : print 'after check_integrity'
+
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.75, y0=0.0, y1=1.0, h=5.0, h0=0.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.y0 = y0
+        self.y1 = y1
+        self.h  = h
+        self.h0 = h0
+
+    def __call__(self, x, y):
+        return self.h0 + self.h*((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1))
+
+
+
+#domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
+
+Set_quantity(domain,'stage',center=(0.0,0.0), radius=0.5, value = 2.0)
+
+if myid == 0 : print 'after set quantity'
+
+
+
+yieldstep = 0.005
+finaltime = 0.05
+
+barrier()
+
+t0 = time.time()
+
+#===========================================================================
+# Main Evolve Loop
+#===========================================================================
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+        sys.stdout.flush()
+        
+        
+
+
+for p in range(numprocs):
+    barrier()
+    if myid == p:
+        print 50*'='
+        print 'P%g' %(myid)
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'Communication time %.2f seconds'%domain.communication_time
+        print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+        print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+        sys.stdout.flush()
+
+
+if domain.number_of_global_triangles < 50000:
+    if myid == 0 :
+        print 'Create dump of triangulation for %g triangles' % domain.number_of_global_triangles
+    domain.dump_triangulation(filename="rectangular_cross_%g.png"% numprocs)
+
+finalize()
Index: /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangular_cross_profile.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangular_cross_profile.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangular_cross_profile.py	(revision 9507)
@@ -0,0 +1,182 @@
+
+"""Run parallel shallow water domain.
+
+   run using command like:
+
+   mpirun -np m python run_parallel_sw_rectangular_cross.py
+
+   where m is the number of processors to be used.
+   
+   Will produce sww files with names domain_Pn_m.sww where n is number of processors and
+   m in [0, n-1] refers to specific processor that owned this part of the partitioned mesh.
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import os
+import sys
+import time
+import numpy as num
+
+#------------------------
+# ANUGA Modules
+#------------------------
+	
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+from anuga import Transmissive_boundary
+
+from anuga import rectangular_cross_domain
+from anuga import create_domain_from_file
+
+
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
+
+
+#--------------------------------------------------------------------------
+# Setup parameters
+#--------------------------------------------------------------------------
+
+#mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0
+#mesh_filename = "merimbula_43200.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
+#mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
+
+
+verbose = True
+
+#--------------------------------------------------------------------------
+# Setup procedures
+#--------------------------------------------------------------------------
+class Set_Stage:
+    """Set an initial condition with constant water height, for x0<x<x1
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+#--------------------------------------------------------------------------
+# Setup Domain only on processor 0
+#--------------------------------------------------------------------------
+if myid == 0:
+    length = 2.0
+    width = 2.0
+    dx = dy = 0.01
+    domain = rectangular_cross_domain(int(length/dx), int(width/dy),
+                                              len1=length, len2=width)
+
+    domain.set_quantity('elevation', -1.0)
+    domain.set_quantity('stage', Set_Stage(h=2.0))
+else:
+    domain = None
+
+#--------------------------------------------------------------------------
+# Distribute sequential domain on processor 0 to other processors
+#--------------------------------------------------------------------------
+
+if myid == 0 and verbose: print 'DISTRIBUTING DOMAIN'
+domain = distribute(domain)
+
+#domain.smooth = False
+barrier()
+for p in range(numprocs):
+    if myid == p:
+        print 'Process ID %g' %myid
+        print 'Number of triangles %g ' %domain.get_number_of_triangles()
+        sys.stdout.flush()
+        
+    barrier()
+
+
+domain.set_flow_algorithm(2.0)
+
+if myid == 0:
+    domain.print_algorithm_parameters()
+    sys.stdout.flush()
+    
+barrier()
+
+domain.set_name('rectangular_cross')
+domain.set_store(False)
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+# This must currently happen *after* domain has been distributed
+#------------------------------------------------------------------------------
+Br = Reflective_boundary(domain)      # Solid reflective wall
+
+domain.set_boundary({'top' :Br, 'bottom' :Br, 'left' :Br, 'right' :Br, 'ghost' : None })
+
+"""
+barrier()
+for p in range(numprocs):
+    if myid == p:
+        print domain.boundary_statistics()
+        sys.stdout.flush()
+        
+    barrier()
+"""
+
+
+domain.dump_triangulation()
+
+#------------------------------------------------------------------------------
+# Evolution
+#------------------------------------------------------------------------------
+if myid == 0 and verbose: print 'EVOLVE'
+
+t0 = time.time()
+finaltime = 0.25
+yieldstep = 0.05
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+
+
+
+## Profiling
+#import cProfile
+#prof_file = 'evolve-prof'+ str(numprocs) + '_' + str(myid) + '.dat'
+#cProfile.run(s,prof_file)
+
+
+barrier()
+
+#for id in range(numprocs):
+#    if myid == id:
+#        import pstats
+#        p = pstats.Stats(prof_file)
+#        #p.sort_stats('cumulative').print_stats(25)
+#        p.sort_stats('time').print_stats(25)
+#        sys.stdout.flush
+#
+#    barrier()
+#
+#barrier()
+
+for p in range(numprocs):
+    if myid == p:
+        print 'Process ID %g' %myid
+        print 'Number of processors %g ' %numprocs
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'Communication time %.2f seconds'%domain.communication_time
+        print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+        print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+
+    barrier()
+
+
+finalize()
+
+
+
+
Index: /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangular_tsunami.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangular_tsunami.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_parallel_sw_rectangular_tsunami.py	(revision 9507)
@@ -0,0 +1,159 @@
+#########################################################
+#
+#  Main file for parallel mesh testing.
+#
+#  This is a modification of the run_parallel_advection.py
+# file.
+#
+#
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+
+import time
+import sys
+import numpy
+import anuga
+
+
+#----------------------------
+# Parallel interface
+#---------------------------
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
+
+
+#--------------------------------------------------------------------------
+# Setup functions for topograpy etc
+#--------------------------------------------------------------------------
+scale_me=1.0
+
+def topography(x,y):
+	return (-x/2.0 +0.05*numpy.sin((x+y)*200.0))*scale_me
+
+def stagefun(x,y):
+    stge=-0.2*scale_me #+0.01*(x>0.9)
+    #topo=topography(x,y)
+    return stge#*(stge>topo) + (topo)*(stge<=topo)
+
+
+
+#--------------------------------------------------------------------------
+# Create domains
+#--------------------------------------------------------------------------
+
+
+
+t0 = time.time()
+
+verbose = True
+
+#--------------------------------------------------------------------------
+# Setup Domain only on processor 0
+#--------------------------------------------------------------------------
+if myid == 0:
+    length = 2.0
+    width = 2.0
+    dx = dy = 0.005
+    #dx = dy = 0.00125
+    dx = dy  = 0.05
+    domain = anuga.rectangular_cross_domain(int(length/dx), int(width/dy),
+                                              len1=length, len2=width, verbose=verbose)
+
+    #---------------------------------------
+    # Add these two commands to use Gareth's
+    # tsunami algorithm. Play with the 
+    # minimum allowed height to remove possible 
+    # unrealistic large velocities
+    #---------------------------------------
+    domain.set_flow_algorithm('tsunami')
+    domain.set_minimum_allowed_height(0.01)
+
+    domain.set_store(True)
+    domain.set_quantity('elevation',topography)     # Use function for elevation
+    domain.get_quantity('elevation').smooth_vertex_values()
+    domain.set_quantity('friction',0.03)            # Constant friction
+    domain.set_quantity('stage', stagefun)          # Constant negative initial stage
+    domain.get_quantity('stage').smooth_vertex_values()
+
+    domain.set_name('rectangular_tsunami')
+
+    domain.print_statistics()
+else:
+    domain = None
+
+t1 = time.time()
+
+if myid == 0 :
+    print 'Create sequential domain ',t1-t0
+
+if myid == 0 and verbose: 
+    print 'DISTRIBUTING DOMAIN'
+    sys.stdout.flush()
+    
+barrier()
+
+#-------------------------------------------------------------------------
+# Distribute domain
+#-------------------------------------------------------------------------
+domain = distribute(domain,verbose=verbose)
+
+
+t2 = time.time()
+
+if myid == 0 :
+    print 'Distribute domain ',t2-t1
+    
+if myid == 0 : print 'after parallel domain'
+
+
+#Boundaries
+T = anuga.Transmissive_boundary(domain)
+R = anuga.Reflective_boundary(domain)
+D = anuga.Dirichlet_boundary([-0.1*scale_me,0.,0.])
+
+
+domain.set_boundary( {'left': R, 'right': D, 'bottom': R, 'top': R, 'ghost': None} )
+
+
+if myid == 0 : print 'after set_boundary'
+
+
+yieldstep = 0.2
+finaltime = 20.0
+
+barrier()
+
+t0 = time.time()
+
+#===========================================================================
+# Main Evolve Loop
+#===========================================================================
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+        sys.stdout.flush()
+						
+
+
+
+for p in range(numprocs):
+    barrier()
+    if myid == p:
+        print 50*'='
+        print 'P%g' %(myid)
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'Communication time %.2f seconds'%domain.communication_time
+        print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+        print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+        sys.stdout.flush()
+
+
+if domain.number_of_global_triangles < 50000:
+    if myid == 0 :
+        print 'Create dump of triangulation for %g triangles' % domain.number_of_global_triangles
+    domain.dump_triangulation(filename="rectangular_cross_%g.png"% numprocs)
+
+finalize()
Index: /trunk/anuga_core/demos/parallel/run_sequential_dist_all_merimbula.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_sequential_dist_all_merimbula.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_sequential_dist_all_merimbula.py	(revision 9507)
@@ -0,0 +1,159 @@
+"""Run parallel shallow water domain.
+
+   run using command like:
+
+   mpirun -np m python run_sequential_dist_all.py
+
+   where m is the number of processors to be used.
+   
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import os
+import sys
+import time
+import numpy as num
+
+#------------------------
+# ANUGA Modules
+#------------------------
+	
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+from anuga import Transmissive_boundary
+
+from anuga import rectangular_cross
+from anuga import create_domain_from_file
+
+
+#from anuga_parallel.sequential_distribute import sequential_distribute
+
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
+
+
+#--------------------------------------------------------------------------
+# Setup parameters
+#--------------------------------------------------------------------------
+
+#mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0
+mesh_filename = "merimbula_17156.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
+#mesh_filename = "merimbula_43200_1.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
+#mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
+#mesh_filename = "test-20.tsh" ; x0 = 250.0 ; x1 = 350.0
+yieldstep = 50
+finaltime = 1500
+verbose = True
+
+#--------------------------------------------------------------------------
+# Setup procedures
+#--------------------------------------------------------------------------
+class Set_Stage:
+    """Set an initial condition with constant water height, for x0<x<x1
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))+1.0
+
+
+class Set_Elevation:
+    """Set an elevation
+    """
+
+    def __init__(self, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return x/self.h
+    
+
+#--------------------------------------------------------------------------
+# Setup Sequential Domain
+#--------------------------------------------------------------------------
+if myid == 0 :
+    domain = create_domain_from_file(mesh_filename)
+    domain.set_quantity('stage', Set_Stage(x0, x1, 2.0))
+    #domain.set_datadir('.')
+    domain.set_name('merimbula_new')
+    domain.set_store(True)
+else:
+    domain = None
+
+
+#--------------------------------------------------------------------------
+# Distribute sequential domain on processor 0 to other processors
+#--------------------------------------------------------------------------
+if verbose: print 'DISTRIBUTING DOMAIN'
+domain = distribute(domain)
+
+
+#--------------------------------------------------------------------------
+# On all processors, setup evolve parameters for domains on all processors
+# (all called "domain"
+#--------------------------------------------------------------------------
+
+domain.set_flow_algorithm('DE0')
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+# This must currently happen *after* domain has been distributed
+#------------------------------------------------------------------------------
+Br = Reflective_boundary(domain)      # Solid reflective wall
+
+domain.set_boundary({'outflow' :Br, 'inflow' :Br, 'inner' :Br, 'exterior' :Br, 'open' :Br})
+
+#------------------------------------------------------------------------------
+# Evolution
+#------------------------------------------------------------------------------
+if myid == 0 and verbose: print 'EVOLVE'
+
+barrier()
+t0 = time.time()
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+
+
+barrier()
+
+for p in range(numprocs):
+    if myid == p:
+        print 'Processor %g ' %myid
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'Communication time %.2f seconds'%domain.communication_time
+        print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+        print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+    else:
+        pass
+
+    barrier()
+
+
+#--------------------------------------------------
+# Merge the individual sww files into one file
+#--------------------------------------------------
+domain.sww_merge(delete_old=False)
+
+finalize()
+
+
+
+
+
+
+
+
+
+
+
Index: /trunk/anuga_core/demos/parallel/run_sequential_dist_distribute_merimbula.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_sequential_dist_distribute_merimbula.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_sequential_dist_distribute_merimbula.py	(revision 9507)
@@ -0,0 +1,100 @@
+
+"""Run parallel shallow water domain.
+
+   run using command like:
+
+   python run_sequential_dist_distribute_merimbula.py
+
+   will produce 4 pickled files corresponding partitioned domains
+   
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import os
+import sys
+import time
+import numpy as num
+
+#------------------------
+# ANUGA Modules
+#------------------------
+	
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+from anuga import Transmissive_boundary
+
+from anuga import rectangular_cross
+from anuga import create_domain_from_file
+
+
+from anuga_parallel.sequential_distribute import sequential_distribute_dump
+
+
+#--------------------------------------------------------------------------
+# Setup parameters
+#--------------------------------------------------------------------------
+
+#mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0
+mesh_filename = "merimbula_17156.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
+#mesh_filename = "merimbula_43200_1.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
+#mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
+#mesh_filename = "test-20.tsh" ; x0 = 250.0 ; x1 = 350.0
+yieldstep = 50
+finaltime = 1500
+verbose = True
+
+#--------------------------------------------------------------------------
+# Setup procedures
+#--------------------------------------------------------------------------
+class Set_Stage:
+    """Set an initial condition with constant water height, for x0<x<x1
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))+1.0
+
+
+class Set_Elevation:
+    """Set an elevation
+    """
+
+    def __init__(self, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return x/self.h
+    
+
+#--------------------------------------------------------------------------
+# Setup Sequential Domain
+#--------------------------------------------------------------------------
+domain = create_domain_from_file(mesh_filename)
+domain.set_quantity('stage', Set_Stage(x0, x1, 2.0))
+#domain.set_datadir('.')
+domain.set_name('merimbula_new')
+domain.set_store(True)
+
+
+#--------------------------------------------------------------------------
+# Distribute sequential domain on processor 0 to other processors
+#--------------------------------------------------------------------------
+
+if verbose: print 'DISTRIBUTING DOMAIN'
+sequential_distribute_dump(domain, 4, verbose=True)
+
+
+
+
+
Index: /trunk/anuga_core/demos/parallel/run_sequential_dist_distribute_rectangular.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_sequential_dist_distribute_rectangular.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_sequential_dist_distribute_rectangular.py	(revision 9507)
@@ -0,0 +1,96 @@
+#########################################################
+#
+#  Main file for parallel mesh testing.
+#
+#  This is a modification of the run_parallel_advection.py
+# file.
+#
+#
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+
+import time
+import sys
+import numpy
+
+
+#----------------------------
+# Sequential interface
+#---------------------------
+from anuga import Transmissive_boundary, Reflective_boundary
+from anuga import rectangular_cross_domain
+
+#----------------------------
+# Parallel interface
+#---------------------------
+from anuga_parallel.sequential_distribute import sequential_distribute_dump
+
+
+t0 = time.time()
+
+verbose = True
+
+#--------------------------------------------------------------------------
+# Setup functions for topograpy etc
+#--------------------------------------------------------------------------
+scale_me=1.0
+
+def topography(x,y):
+	return (-x/2.0 +0.05*numpy.sin((x+y)*200.0))*scale_me
+
+def stagefun(x,y):
+    stge=-0.2*scale_me #+0.01*(x>0.9)
+    #topo=topography(x,y)
+    return stge#*(stge>topo) + (topo)*(stge<=topo)
+
+
+#--------------------------------------------------------------------------
+# Setup Domain only on processor 0
+#--------------------------------------------------------------------------
+myid = 0
+numprocs = 3
+
+
+
+length = 2.0
+width = 2.0
+dx = dy = 0.005  # 640,000
+dx = dy  = 0.05
+domain = rectangular_cross_domain(int(length/dx), int(width/dy),
+                                          len1=length, len2=width, verbose=verbose)
+
+
+print domain.number_of_global_triangles
+
+domain.set_store(True)
+domain.set_flow_algorithm('DE0')
+domain.set_minimum_allowed_height(0.01)
+domain.set_quantity('elevation',topography)     # Use function for elevation
+domain.get_quantity('elevation').smooth_vertex_values()
+domain.set_quantity('friction',0.03)            # Constant friction
+domain.set_quantity('stage', stagefun)          # Constant negative initial stage
+domain.get_quantity('stage').smooth_vertex_values()
+
+domain.set_name('sw_rectangle')
+
+
+t1 = time.time()
+print 'Create sequential domain: time ',t1-t0
+    
+
+#-------------------------------------------------------------------------
+# Distribute domain
+#----------------------------------------------------------------------
+sequential_distribute_dump(domain, numprocs, verbose = True)
+
+
+t2 = time.time()
+print 'Sequential Distribute domain: time ',t2-t1
+    
+
+
+
Index: /trunk/anuga_core/demos/parallel/run_sequential_dist_evolve_merimbula.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_sequential_dist_evolve_merimbula.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_sequential_dist_evolve_merimbula.py	(revision 9507)
@@ -0,0 +1,109 @@
+
+"""Run parallel shallow water domain.
+
+   run using command like:
+
+   mpirun -np 4  python run_sequential_dist_evolve_merimbula.py
+
+   Need to have run run_sequential_dist_distribute_merimbula.py
+
+   first to produce the required pickle files encoding the partitioning
+   """
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import os
+import sys
+import time
+import numpy as num
+
+#------------------------
+# ANUGA Modules
+#------------------------
+	
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Dirichlet_boundary
+from anuga import Time_boundary
+from anuga import Transmissive_boundary
+
+from anuga import rectangular_cross
+from anuga import create_domain_from_file
+
+
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
+
+from anuga_parallel.sequential_distribute import sequential_distribute_load
+
+#--------------------------------------------------------------------------
+# Setup parameters
+#--------------------------------------------------------------------------
+yieldstep = 50
+finaltime = 1500
+verbose = True
+
+
+
+#-------------------------------------------------------------------------------
+# Read in cPickled distribute mesh data
+#-------------------------------------------------------------------------------
+domain = sequential_distribute_load(filename='merimbula_new', verbose = True)
+
+#--------------------------------------------------------------------------
+# On all processors, setup evolve parameters for domains on all processors
+# (all called "domain"
+#--------------------------------------------------------------------------
+
+domain.set_flow_algorithm('DE0')
+
+
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+# This must currently happen *after* domain has been distributed
+#------------------------------------------------------------------------------
+Br = Reflective_boundary(domain)      # Solid reflective wall
+
+domain.set_boundary({'outflow' :Br, 'inflow' :Br, 'inner' :Br, 'exterior' :Br, 'open' :Br})
+
+#------------------------------------------------------------------------------
+# Evolution
+#------------------------------------------------------------------------------
+if myid == 0 and verbose: print 'EVOLVE'
+
+barrier()
+t0 = time.time()
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+
+
+barrier()
+
+for p in range(numprocs):
+    if myid == p:
+        print 'Processor %g ' %myid
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'Communication time %.2f seconds'%domain.communication_time
+        print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+        print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+    else:
+        pass
+
+    barrier()
+
+
+#--------------------------------------------------
+# Merge the individual sww files into one file
+#--------------------------------------------------
+
+
+
+
+finalize()
+
+
+
+
Index: /trunk/anuga_core/demos/parallel/run_sequential_dist_parallel_evolve_rectangular.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_sequential_dist_parallel_evolve_rectangular.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_sequential_dist_parallel_evolve_rectangular.py	(revision 9507)
@@ -0,0 +1,134 @@
+#########################################################
+#
+#  Main file for parallel mesh testing.
+#
+#  This is a modification of the run_parallel_advection.py
+# file.
+#
+#
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+
+import time
+import sys
+import numpy
+import anuga
+
+
+#----------------------------
+# Parallel interface
+#---------------------------
+from anuga_parallel import myid, numprocs, finalize, barrier
+
+from anuga_parallel.sequential_distribute import sequential_distribute_load
+
+
+#--------------------------------------------------------------------------
+# Setup functions for topograpy etc
+#--------------------------------------------------------------------------
+scale_me=1.0
+
+def topography(x,y):
+	return (-x/2.0 +0.05*numpy.sin((x+y)*200.0))*scale_me
+
+def stagefun(x,y):
+    stge=-0.2*scale_me #+0.01*(x>0.9)
+    #topo=topography(x,y)
+    return stge#*(stge>topo) + (topo)*(stge<=topo)
+
+
+
+
+verbose = True
+basename = 'sw_rectangle'
+#basename = 'test'
+
+#--------------------------------------------------------------------------
+# Create domains
+#--------------------------------------------------------------------------
+#import cPickle
+#pickle_name = pickle_base+'_P%g_%g.pickle'% (numprocs,myid)
+#f = file(pickle_name, 'rb')
+# = cPickle.load(f)
+#f.close()
+
+
+t1 = time.time()
+
+if myid == 0 and verbose:
+    print 'Loading DISTRIBUTING DOMAIN'
+    sys.stdout.flush()
+
+#-------------------------------------------------------------------------
+# Distribute domain
+#-------------------------------------------------------------------------
+domain = sequential_distribute_load(filename = basename, verbose=verbose)
+
+
+t2 = time.time()
+
+if myid == 0 :
+    print 'Load Distribute domain ',t2-t1
+    
+if myid == 0 : print 'after parallel domain'
+
+
+#Boundaries
+T = anuga.Transmissive_boundary(domain)
+R = anuga.Reflective_boundary(domain)
+D = anuga.Dirichlet_boundary([-0.1*scale_me,0.,0.])
+
+
+domain.set_boundary( {'left': R, 'right': D, 'bottom': R, 'top': R, 'ghost': None} )
+
+
+if myid == 0 : print 'after set_boundary'
+
+
+yieldstep = 0.2
+finaltime = 20.0
+
+barrier()
+
+t0 = time.time()
+
+#===========================================================================
+# Main Evolve Loop
+#===========================================================================
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    if myid == 0:
+        domain.write_time()
+        sys.stdout.flush()
+						
+
+
+
+for p in range(numprocs):
+    barrier()
+    if myid == p:
+        print 50*'='
+        print 'P%g' %(myid)
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'Communication time %.2f seconds'%domain.communication_time
+        print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time
+        print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time
+        sys.stdout.flush()
+
+
+if domain.number_of_global_triangles < 50000:
+    if myid == 0 :
+        print 'Create dump of triangulation for %g triangles' % domain.number_of_global_triangles
+    domain.dump_triangulation(filename="rectangular_cross_%g.png"% numprocs)
+
+
+#--------------------------------------------------
+# Merge the individual sww files into one file
+#--------------------------------------------------
+domain.sww_merge(delete_old=False)
+
+
+finalize()
Index: /trunk/anuga_core/demos/parallel/run_sw_lost_mass.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_sw_lost_mass.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_sw_lost_mass.py	(revision 9507)
@@ -0,0 +1,108 @@
+#!/usr/bin/env/python
+#########################################################
+#
+# Shows problem of conservation of mass. Due to material being
+# removed if too shallow. Can lead to problems if timestep
+# set very small.
+#
+#
+#
+#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
+# June 2005
+#
+#
+#
+#########################################################
+
+
+##############################################
+# Set a small yieldstep (say 0.005) to see
+# lose of conservation
+# Set it to a larger value (say 0.1) and have
+# conservation "restored"
+##############################################
+yieldstep = 0.005
+#yieldstep = 0.08
+#yieldstep = 0.05
+finaltime = 100.0
+
+import sys
+import time
+
+
+from Numeric import array, zeros, Float
+
+#from shallow_water import Domain
+
+from anuga.shallow_water import Domain
+
+# mesh partition routines
+
+from anuga.abstract_2d_finite_volumes.pmesh2domain\
+     import pmesh_to_domain_instance
+
+# read in the processor information
+
+
+#-------
+# Domain
+rect = zeros( 4, Float) # Buffer for results
+
+class Set_Stage:
+    """Set an initial condition with constant water height, for x0<x<x1
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+
+
+    # read in the test files
+
+filename = 'test-100.tsh'
+
+nx = 1
+ny = 1
+
+domain = pmesh_to_domain_instance(filename, Domain)
+
+domain.set_quantity('stage', Set_Stage(200.0,300.0,5.0))
+rect = array(domain.xy_extent, Float)
+
+try:
+    domain.initialise_visualiser(rect=rect)
+    #domain.visualiser.coloring['stage'] = True
+    domain.visualiser.scale_z['stage'] = 0.2
+    domain.visualiser.scale_z['elevation'] = 0.05
+except:
+    print 'No visualiser'
+
+
+domain.default_order = 2
+
+#Boundaries
+from shallow_water import Transmissive_boundary, Reflective_boundary
+
+T = Transmissive_boundary(domain)
+R = Reflective_boundary(domain)
+domain.set_boundary( {'outflow': R, 'inflow': R, 'inner':R, 'exterior': R, 'open':R} )
+
+
+
+#---------
+# Evolution
+t0 = time.time()
+
+print 'No of triangles %d'%(domain.number_of_triangles)
+
+
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    #domain.write_time()
+    print '    Integral stage = ', domain.quantities['stage'].get_integral(),' Time = ',domain.time
+
+
Index: /trunk/anuga_core/demos/parallel/run_sw_lost_mass_simpler.py
===================================================================
--- /trunk/anuga_core/demos/parallel/run_sw_lost_mass_simpler.py	(revision 9507)
+++ /trunk/anuga_core/demos/parallel/run_sw_lost_mass_simpler.py	(revision 9507)
@@ -0,0 +1,74 @@
+#!/usr/bin/env python
+#########################################################
+#
+# Shows problem of conservation of mass. Due to material being
+# removed if too shallow. Can lead to problems if timestep
+# set very small.
+#
+#
+#
+#  Authors: Linda Stals, Steve Roberts, Matthew Hardy, Ole Nielsen
+#  July 2005
+#
+#
+#########################################################
+
+
+##############################################
+# Change min_depth and see how larger values aggravate the problem
+yieldstep = 0.1
+finaltime = 10.0
+#min_depth = 1.0e-4
+min_depth = 1.0e-2
+
+
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
+from anuga.shallow_water import Domain, Reflective_boundary
+
+
+
+#Create shallow water domain
+points, vertices, boundary = rectangular(50, 50, len1=500, len2=500)
+domain = Domain(points, vertices, boundary)
+domain.smooth = False
+domain.visualise = True
+domain.default_order = 1
+domain.minimum_allowed_height = min_depth
+print 'Extent', domain.get_extent()
+
+# Set initial condition
+class Set_IC:
+    """Set an initial condition with a constant value, for x0<x<x1
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+
+domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
+
+try:
+    domain.initialise_visualiser()
+except:
+    print 'No visualiser'
+else:    
+    domain.visualiser.scale_z['stage'] = 0.2
+    domain.visualiser.scale_z['elevation'] = 0.05
+    
+
+#Boundaries
+R = Reflective_boundary(domain)
+domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R} )
+
+
+# Evolution
+print 'Minimal allowed water height = ', domain.minimum_allowed_height 
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    print 'Integral stage = ', domain.quantities['stage'].get_integral(),' Time = ',domain.time
+
+
Index: /trunk/anuga_core/demos/structures/run_culvert.py
===================================================================
--- /trunk/anuga_core/demos/structures/run_culvert.py	(revision 9507)
+++ /trunk/anuga_core/demos/structures/run_culvert.py	(revision 9507)
@@ -0,0 +1,162 @@
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.geometry.polygon_function import Polygon_function
+        
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+import anuga
+
+from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
+from anuga.structures.inlet_operator import Inlet_operator
+                            
+#from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+     
+from math import pi, pow, sqrt
+
+import numpy as num
+
+
+
+
+
+"""test_that_culvert_runs_rating
+
+This test exercises the culvert and checks values outside rating curve
+are dealt with       
+"""
+
+path = get_pathname_from_package('anuga.culvert_flows')    
+
+length = 40.
+width = 15.
+
+dx = dy = 0.5          # Resolution: Length of subdivisions on both axes
+
+points, vertices, boundary = rectangular_cross(int(length/dx),
+                                               int(width/dy),
+                                               len1=length, 
+                                               len2=width)
+domain = anuga.Domain(points, vertices, boundary)   
+domain.set_starttime(10)
+domain.set_name()                 # Output name
+domain.set_default_order(2)
+#domain.set_beta(1.5)
+
+
+#----------------------------------------------------------------------
+# Setup initial conditions
+#----------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=-x/1000
+    
+    N = len(x)
+    for i in range(N):
+
+       # Sloping Embankment Across Channel
+        if 5.0 < x[i] < 10.1:
+            # Cut Out Segment for Culvert face                
+            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+               z[i]=z[i]
+            else:
+               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+        if 10.0 < x[i] < 12.1:
+           z[i] +=  2.5                    # Flat Crest of Embankment
+        if 12.0 < x[i] < 14.5:
+            # Cut Out Segment for Culvert face                
+            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+               z[i]=z[i]
+            else:
+               z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                   
+        
+    return z
+
+
+domain.set_quantity('elevation', topography) 
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+filename=os.path.join(path, 'example_rating_curve.csv')
+
+end_point0 = num.array([9.0, 2.5])
+end_point1 = num.array([13.0, 2.5])
+
+Boyd_pipe_operator(domain,
+                    #end_point0=[9.0, 2.5],
+                    #end_point1=[13.0, 2.5],
+                    #exchange_line0=[[9.0, 1.75],[9.0, 3.25]],
+                    #exchange_line1=[[13.0, 1.75],[13.0, 3.25]],
+                    losses=1.5,
+                    end_points=[end_point0, end_point1],
+                    diameter=1.5,
+                    apron=0.5,
+                    use_momentum_jet=True, 
+                    use_velocity_head=False,
+                    manning=0.013,
+                    verbose=False)
+
+
+line = [[0.0, 5.0], [0.0, 10.0]]
+Q = 5.0
+#Inlet_operator(domain, line, Q)
+
+
+
+
+##-----------------------------------------------------------------------
+## Setup boundary conditions
+##-----------------------------------------------------------------------
+
+## Inflow based on Flow Depth and Approaching Momentum
+Bi = anuga.Dirichlet_boundary([2.0, 0.0, 0.0])
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+#Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+## Upstream and downstream conditions that will exceed the rating curve
+## I.e produce delta_h outside the range [0, 10] specified in the the 
+## file example_rating_curve.csv
+#Btus = anuga.Time_boundary(domain, \
+            #lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
+#Btds = anuga.Time_boundary(domain, \
+            #lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
+#domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+##-----------------------------------------------------------------------
+## Evolve system through time
+##-----------------------------------------------------------------------
+
+#min_delta_w = sys.maxint 
+#max_delta_w = -min_delta_w
+for t in domain.evolve(yieldstep=1.0, finaltime=50.0):
+    domain.write_time()
+
+    #if domain.get_time() > 150.5 and domain.ge t_time() < 151.5 :
+        #Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
+        #domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
+
+    #delta_w = culvert.inlet.stage - culvert.outlet.stage
+    
+    #if delta_w > max_delta_w: max_delta_w = delta_w
+    #if delta_w < min_delta_w: min_delta_w = delta_w            
+    
+    pass
+
+## Check that extreme values in rating curve have been exceeded
+## so that we know that condition has been exercised
+#assert min_delta_w < 0
+#assert max_delta_w > 10        
+
+
+#os.remove('Test_culvert.sww')
Index: /trunk/anuga_core/demos/structures/run_culvert_inlet.py
===================================================================
--- /trunk/anuga_core/demos/structures/run_culvert_inlet.py	(revision 9507)
+++ /trunk/anuga_core/demos/structures/run_culvert_inlet.py	(revision 9507)
@@ -0,0 +1,158 @@
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.geometry.polygon_function import Polygon_function
+        
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+import anuga
+
+from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
+from anuga.structures.inlet_operator import Inlet_operator
+                            
+#from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+     
+from math import pi, pow, sqrt
+
+import numpy as num
+
+
+
+
+
+"""test_that_culvert_runs_rating
+
+This test exercises the culvert and checks values outside rating curve
+are dealt with       
+"""
+
+path = get_pathname_from_package('anuga.culvert_flows')    
+
+length = 40.
+width = 15.
+
+dx = dy = 0.5          # Resolution: Length of subdivisions on both axes
+
+points, vertices, boundary = rectangular_cross(int(length/dx),
+                                               int(width/dy),
+                                               len1=length, 
+                                               len2=width)
+domain = anuga.Domain(points, vertices, boundary)   
+domain.set_name('run_culvert_inlet')                 # Output name
+domain.set_default_order(2)
+#domain.set_beta(1.5)
+
+
+#----------------------------------------------------------------------
+# Setup initial conditions
+#----------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=-x/1000
+    
+    N = len(x)
+    for i in range(N):
+
+       # Sloping Embankment Across Channel
+        if 5.0 < x[i] < 10.1:
+            # Cut Out Segment for Culvert face                
+            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+               z[i]=z[i]
+            else:
+               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+        if 10.0 < x[i] < 12.1:
+           z[i] +=  2.5                    # Flat Crest of Embankment
+        if 12.0 < x[i] < 14.5:
+            # Cut Out Segment for Culvert face                
+            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+               z[i]=z[i]
+            else:
+               z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                   
+        
+    return z
+
+
+domain.set_quantity('elevation', topography) 
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+filename=os.path.join(path, 'example_rating_curve.csv')
+
+
+
+Boyd_pipe_operator(domain,
+                            end_points=[[9.0, 2.5],[13.0, 2.5]],
+                            losses=1.5,
+                            diameter=1.5,
+                            apron=5.0,
+                            use_momentum_jet=True,
+                            use_velocity_head=False,
+                            manning=0.013,
+                            verbose=False)
+
+
+line = [[0.0, 5.0], [0.0, 10.0]]
+Q = 5.0
+Inlet_operator(domain, line, Q)
+
+
+
+
+##-----------------------------------------------------------------------
+## Setup boundary conditions
+##-----------------------------------------------------------------------
+
+## Inflow based on Flow Depth and Approaching Momentum
+Bi = anuga.Dirichlet_boundary([2.0, 0.0, 0.0])
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+#Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+## Upstream and downstream conditions that will exceed the rating curve
+## I.e produce delta_h outside the range [0, 10] specified in the the 
+## file example_rating_curve.csv
+#Btus = anuga.Time_boundary(domain, \
+            #lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
+#Btds = anuga.Time_boundary(domain, \
+            #lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
+#domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+##-----------------------------------------------------------------------
+## Evolve system through time
+##-----------------------------------------------------------------------
+
+#min_delta_w = sys.maxint 
+#max_delta_w = -min_delta_w
+for t in domain.evolve(yieldstep = 1.0, finaltime = 200):
+    domain.write_time()
+
+    #if domain.get_time() > 150.5 and domain.get_time() < 151.5 :
+        #Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
+        #domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
+
+    #delta_w = culvert.inlet.stage - culvert.outlet.stage
+    
+    #if delta_w > max_delta_w: max_delta_w = delta_w
+    #if delta_w < min_delta_w: min_delta_w = delta_w
+
+    print domain.volumetric_balance_statistics()
+    
+    pass
+
+## Check that extreme values in rating curve have been exceeded
+## so that we know that condition has been exercised
+#assert min_delta_w < 0
+#assert max_delta_w > 10        
+
+
+#os.remove('Test_culvert.sww')
Index: /trunk/anuga_core/demos/structures/run_gate_operator.py
===================================================================
--- /trunk/anuga_core/demos/structures/run_gate_operator.py	(revision 9507)
+++ /trunk/anuga_core/demos/structures/run_gate_operator.py	(revision 9507)
@@ -0,0 +1,146 @@
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.geometry.polygon_function import Polygon_function
+        
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+import anuga
+
+from anuga.structures.boyd_box_operator import Boyd_box_operator
+from anuga.structures.inlet_operator import Inlet_operator
+                            
+#from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+     
+from math import pi, pow, sqrt
+
+import numpy as num
+
+
+
+
+
+"""test_that_culvert_runs_rating
+
+This test exercises the culvert and checks values outside rating curve
+are dealt with       
+"""
+
+path = get_pathname_from_package('anuga.culvert_flows')    
+
+length = 40.
+width = 15.
+
+dx = dy = 0.5          # Resolution: Length of subdivisions on both axes
+
+points, vertices, boundary = rectangular_cross(int(length/dx),
+                                               int(width/dy),
+                                               len1=length, 
+                                               len2=width)
+domain = anuga.Domain(points, vertices, boundary)   
+domain.set_name('run_gate_operator')                 # Output name
+domain.set_default_order(2)
+#domain.set_beta(1.5)
+
+
+#----------------------------------------------------------------------
+# Setup initial conditions
+#----------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=-x/1000
+    
+    N = len(x)
+    for i in range(N):
+
+       # Sloping Embankment Across Channel
+        if 5.0 < x[i] < 10.1:
+            # Cut Out Segment for Culvert face                
+            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+               z[i]=z[i]
+            else:
+               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+        if 10.0 < x[i] < 12.1:
+           z[i] +=  2.5                    # Flat Crest of Embankment
+        if 12.0 < x[i] < 14.5:
+            # Cut Out Segment for Culvert face                
+            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+               z[i]=z[i]
+            else:
+               z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                   
+        
+    return z
+
+
+domain.set_quantity('elevation', topography) 
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+filename=os.path.join(path, 'example_rating_curve.csv')
+
+
+gate = Boyd_box_operator(domain,
+                            end_points=[[9.0, 2.5],[13.0, 2.5]],
+                            losses=1.5,
+                            width=1.5,
+                            height = 10.0,
+                            apron=5.0,
+                            use_momentum_jet=True,
+                            use_velocity_head=False,
+                            manning=0.013,
+                            verbose=False)
+
+line = [[0.0, 5.0], [0.0, 10.0]]
+Q = 1.0
+Inlet_operator(domain, line, Q)
+
+
+
+
+##-----------------------------------------------------------------------
+## Setup boundary conditions
+##-----------------------------------------------------------------------
+
+## Inflow based on Flow Depth and Approaching Momentum
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+##-----------------------------------------------------------------------
+## Evolve system through time
+##-----------------------------------------------------------------------
+
+#min_delta_w = sys.maxint 
+#max_delta_w = -min_delta_w
+for t in domain.evolve(yieldstep = 1.0, finaltime = 50):
+    domain.write_time()
+
+    if num.allclose(t, 10.0):
+        gate.set_culvert_height(0.0)
+
+    Q, velocity, depth = gate.discharge_routine()
+
+    print gate.inlets[0].get_enquiry_stage()
+    print gate.inlets[1].get_enquiry_stage()
+    print gate.get_culvert_height()
+
+    print Q
+    print velocity
+    print depth
+        
+    gate.print_timestepping_statistics()
+    
+
+
+
+
Index: /trunk/anuga_core/demos/structures/run_inlet_operator.py
===================================================================
--- /trunk/anuga_core/demos/structures/run_inlet_operator.py	(revision 9507)
+++ /trunk/anuga_core/demos/structures/run_inlet_operator.py	(revision 9507)
@@ -0,0 +1,168 @@
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.geometry.polygon_function import Polygon_function
+        
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+from anuga.abstract_2d_finite_volumes.util import file_function
+
+import anuga
+
+from anuga.structures.boyd_box_operator import Boyd_box_operator
+from anuga.structures.inlet_operator import Inlet_operator
+                            
+#from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+     
+from math import pi, pow, sqrt
+
+import numpy as num
+
+
+
+
+
+"""test_that_culvert_runs_rating
+
+This test exercises the culvert and checks values outside rating curve
+are dealt with       
+"""
+
+path = get_pathname_from_package('anuga.culvert_flows')    
+
+length = 40.
+width = 15.
+
+dx = dy = 0.5          # Resolution: Length of subdivisions on both axes
+
+points, vertices, boundary = rectangular_cross(int(length/dx),
+                                               int(width/dy),
+                                               len1=length, 
+                                               len2=width)
+domain = anuga.Domain(points, vertices, boundary)   
+domain.set_name('Testing_inlet_operator')                 # Output name
+domain.set_default_order(2)
+#domain.set_beta(1.5)
+
+
+#----------------------------------------------------------------------
+# Setup initial conditions
+#----------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=-x/1000
+    
+    N = len(x)
+    for i in range(N):
+
+       # Sloping Embankment Across Channel
+        if 5.0 < x[i] < 10.1:
+            # Cut Out Segment for Culvert face                
+            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+               z[i]=z[i]
+            else:
+               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+        if 10.0 < x[i] < 12.1:
+           z[i] +=  2.5                    # Flat Crest of Embankment
+        if 12.0 < x[i] < 14.5:
+            # Cut Out Segment for Culvert face                
+            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+               z[i]=z[i]
+            else:
+               z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                   
+        
+    return z
+
+
+domain.set_quantity('elevation', topography) 
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+filename=os.path.join(path, 'example_rating_curve.csv')
+
+
+
+#Boyd_box_operator(domain,
+#                            end_points=[[9.0, 2.5],[13.0, 2.5]],
+#                            losses=1.5,
+#                            width=1.5,
+#                            apron=5.0,
+#                            use_momentum_jet=True,
+#                            use_velocity_head=False,
+#                            manning=0.013,
+#                            verbose=False)
+
+
+line0 = [[0.0, 5.0], [0.0, 10.0]]
+
+
+Q = file_function('test_hydrograph.tms', quantities=['hydrograph'])
+
+#inlet0 = Inlet_operator(domain, line0, Q, label='first inlet')
+
+
+line1 = [[1.0, 5.0], [2.0, 10.0]]
+poly1 = [[1.0, 5.0], [2.0, 5.0], [2.0, 10.0], [1.0, 10.0]]
+inlet1 = Inlet_operator(domain, poly1, 2.0)
+
+
+
+##-----------------------------------------------------------------------
+## Setup boundary conditions
+##-----------------------------------------------------------------------
+
+## Inflow based on Flow Depth and Approaching Momentum
+Bi = anuga.Dirichlet_boundary([2.0, 0.0, 0.0])
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+#Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+## Upstream and downstream conditions that will exceed the rating curve
+## I.e produce delta_h outside the range [0, 10] specified in the the 
+## file example_rating_curve.csv
+#Btus = anuga.Time_boundary(domain, \
+            #lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
+#Btds = anuga.Time_boundary(domain, \
+            #lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
+#domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+##-----------------------------------------------------------------------
+## Evolve system through time
+##-----------------------------------------------------------------------
+
+#min_delta_w = sys.maxint 
+#max_delta_w = -min_delta_w
+for t in domain.evolve(yieldstep = 1.0, finaltime = 38):
+    domain.write_time()
+
+    #if domain.get_time() > 150.5 and domain.get_time() < 151.5 :
+        #Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
+        #domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
+
+    #delta_w = culvert.inlet.stage - culvert.outlet.stage
+    
+    #if delta_w > max_delta_w: max_delta_w = delta_w
+    #if delta_w < min_delta_w: min_delta_w = delta_w
+
+    print domain.volumetric_balance_statistics()
+
+    #inlet0.print_timestepping_statistics()
+    inlet1.print_timestepping_statistics()
+    pass
+
+## Check that extreme values in rating curve have been exceeded
+## so that we know that condition has been exercised
+#assert min_delta_w < 0
+#assert max_delta_w > 10        
+
+
+#os.remove('Test_culvert.sww')
Index: /trunk/anuga_core/demos/structures/run_open_slot_wide_bridge.py
===================================================================
--- /trunk/anuga_core/demos/structures/run_open_slot_wide_bridge.py	(revision 9507)
+++ /trunk/anuga_core/demos/structures/run_open_slot_wide_bridge.py	(revision 9507)
@@ -0,0 +1,191 @@
+""" Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction
+
+This example includes a Model Topography that shows a TYPICAL Headwall Configuration
+
+The aim is to change the Culvert Routine to Model more precisely the abstraction
+from a vertical face.
+
+The inflow must include the impact of Approach velocity.
+Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style
+abstraction
+
+"""
+print 'Starting.... Importing Modules...'
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import anuga
+
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water.shallow_water_domain import Domain
+from anuga.shallow_water.forcing import Rainfall, Inflow
+#from anuga.shallow_water.forcing import Reflective_boundary
+#from anuga.shallow_water.forcing import Dirichlet_boundary
+#from anuga.shallow_water.forcing import Transmissive_boundary, Time_boundary
+
+from anuga.culvert_flows.culvert_class import Culvert_flow
+from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+#from anuga.culvert_flows.culvert_routines import weir_orifice_channel_culvert_model
+from math import pi,pow,sqrt
+
+import numpy as num
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+print 'Setting up domain'
+
+length = 200. #x-Dir
+width = 200.  #y-dir
+
+dx = dy = 2.5          # Resolution: Length of subdivisions on both axes
+#dx = dy = .5           # Resolution: Length of subdivisions on both axes
+#dx = dy = .5           # Resolution: Length of subdivisions on both axes
+#dx = dy = .1           # 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)   
+domain.set_name('Test_Open_slot_WIDE_BRIDGE')                 # Output name
+domain.set_default_order(2)
+domain.H0 = 0.01
+domain.tight_slope_limiters = 1
+
+print 'Size', len(domain)
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=10.0-x/100.0  # % Longitudinal Slope
+    
+    #       NOW Add bits and Pieces to topography
+    bank_hgt=10.0
+    bridge_width = 50.0
+    bank_width = 10.0
+    
+    us_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
+    us_start_x = 10.0
+    top_start_y = 50.0
+    us_slope = 3.0  #Horiz : Vertic
+    ds_slope = 3.0
+    ds_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
+    centre_line_y= top_start_y+bridge_width/2.0
+
+    # CALCULATE PARAMETERS TO FORM THE EMBANKMENT
+    us_slope_length = bank_hgt*us_slope
+    us_end_x =us_start_x + us_slope_length
+    us_toe_start_y =top_start_y - us_slope_length / us_apron_skew
+    us_toe_end_y = top_start_y + bridge_width + us_slope_length / us_apron_skew
+
+    top_end_y = top_start_y + bridge_width
+    ds_slope_length = bank_hgt*ds_slope
+    ds_start_x = us_end_x + bank_width
+    ds_end_x = ds_start_x + ds_slope_length
+
+    ds_toe_start_y =top_start_y - ds_slope_length / ds_apron_skew
+    ds_toe_end_y = top_start_y + bridge_width + ds_slope_length / ds_apron_skew
+
+
+    N = len(x)
+    for i in range(N):
+
+       # Sloping Embankment Across Channel
+        if us_start_x < x[i] < us_end_x +0.1:   # For UPSLOPE on the Upstream FACE
+        #if 5.0 < x[i] < 10.1: # For a Range of X, and over a Range of Y based on X adjust Z
+            if us_toe_start_y +(x[i] - us_start_x)/us_apron_skew < y[i] < us_toe_end_y - (x[i] - us_start_x)/us_apron_skew:
+			#if  49.0+(x[i]-5.0)/5.0 <  y[i]  < 151.0 - (x[i]-5.0)/5.0: # Cut Out Base Segment for Culvert FACE
+               z[i]=z[i] # Flat Apron
+            else:
+               z[i] += z[i] + (x[i] - us_start_x)/us_slope    # Set elevation for Sloping Segment  U/S Face
+        if us_end_x < x[i] < ds_start_x + 0.1: # FOR The TOP of BANK Segment
+           if top_start_y < y[i] < top_end_y:
+               z[i]=z[i] # Flat Apron
+           else:
+               z[i] +=  z[i]+bank_hgt        # Flat Crest of Embankment
+        if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face
+            if  top_start_y-(x[i]-ds_start_x)/ds_apron_skew <  y[i]  < top_end_y + (x[i]-ds_start_x)/ds_apron_skew: # Cut Out Segment for Culvert FACE
+               z[i]=z[i] # Flat Apron
+            else:
+               z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope       # Sloping D/S Face
+           
+        
+
+    return z
+
+print 'Setting Quantities....'
+domain.set_quantity('elevation', topography)  # Use function for elevation
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+
+
+
+#------------------------------------------------------------------------------
+# Setup specialised forcing terms
+#------------------------------------------------------------------------------
+
+#------------------------------------------------------------------------------
+# Setup CULVERT INLETS and OUTLETS in Current Topography
+#------------------------------------------------------------------------------
+print 'DEFINING any Structures if Required'
+
+#  DEFINE CULVERT INLET AND OUTLETS
+
+"""
+culvert_rating = Culvert_flow(domain,
+                       culvert_description_filename='example_rating_curve.csv',
+                       end_point0=[40.0, 75.0], 
+                       end_point1=[50.0, 75.0],
+                       verbose=True)
+
+
+culvert_energy = Culvert_flow(domain,
+                       label='Culvert No. 1',
+                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
+                       end_point0=[40.0, 75.0], 
+                       end_point1=[50.0, 75.0],
+                       width=50.0,depth=5.0,
+                       culvert_routine=boyd_generalised_culvert_model,        
+                       number_of_barrels=1,
+                       #update_interval=0.25,
+                       log_file=True,
+                       discharge_hydrograph=True,
+                       use_velocity_head=False,
+                       verbose=True)
+
+domain.forcing_terms.append(culvert_energy)
+"""
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+print 'Setting Boundary Conditions'
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
+
+Bo = anuga.Dirichlet_boundary([-5.0, 0, 0])           # Outflow water at -5.0
+Bd = anuga.Dirichlet_boundary([0,0,0])                # Outflow water at 0.0
+
+#Btus = Time_boundary(domain, lambda t: [0.0+ 1.025*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0])
+#Btds = Time_boundary(domain, lambda t: [0.0+ 0.0075*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0])
+
+Btus = anuga.Dirichlet_boundary([18.0, 0, 0])           # Outflow water at 5.0
+Btds = anuga.Dirichlet_boundary([0.0, 0, 0])           # Outflow water at 1.0
+domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 1, finaltime = 100):
+    print domain.timestepping_statistics()
Index: /trunk/anuga_core/demos/structures/run_outlet_control.py
===================================================================
--- /trunk/anuga_core/demos/structures/run_outlet_control.py	(revision 9507)
+++ /trunk/anuga_core/demos/structures/run_outlet_control.py	(revision 9507)
@@ -0,0 +1,232 @@
+""" Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction
+
+This example includes a Model Topography that shows a TYPICAL Headwall Configuration
+
+The aim is to change the Culvert Routine to Model more precisely the abstraction
+from a vertical face.
+
+The inflow must include the impact of Approach velocity.
+Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style
+abstraction
+
+"""
+print 'Starting.... Importing Modules...'
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import anuga
+
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water.shallow_water_domain import Domain
+from anuga.shallow_water.forcing import Rainfall, Inflow
+#from anuga.shallow_water.forcing import Reflective_boundary
+#from anuga.shallow_water.forcing import Dirichlet_boundary
+#from anuga.shallow_water.forcing import Transmissive_boundary, Time_boundary
+
+#from anuga.culvert_flows.culvert_class import Culvert_flow
+from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
+from anuga.structures.boyd_box_operator import Boyd_box_operator
+#from anuga.culvert_flows.culvert_routines import weir_orifice_channel_culvert_model
+from math import pi,pow,sqrt
+
+import numpy as num
+
+
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+print 'Setting up domain'
+
+length = 120. #x-Dir
+width = 200.  #y-dir
+
+dx = dy = 2.0          # Resolution: Length of subdivisions on both axes
+#dx = dy = .5           # Resolution: Length of subdivisions on both axes
+#dx = dy = .5           # Resolution: Length of subdivisions on both axes
+#dx = dy = .1           # 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)   
+domain.set_name('Test_Outlet_Ctrl')                 # Output name
+domain.set_default_order(2)
+domain.H0 = 0.01
+domain.tight_slope_limiters = 1
+
+print 'Size', len(domain)
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=10.0-x/100.0  # % Longitudinal Slope
+    
+    #       NOW Add bits and Pieces to topography
+    bank_hgt=10.0
+    bridge_width = 50.0
+    bank_width = 10.0
+    
+    us_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
+    us_start_x = 10.0
+    top_start_y = 50.0
+    us_slope = 3.0  #Horiz : Vertic
+    ds_slope = 3.0
+    ds_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
+    centre_line_y= top_start_y+bridge_width/2.0
+
+    # CALCULATE PARAMETERS TO FORM THE EMBANKMENT
+    us_slope_length = bank_hgt*us_slope
+    us_end_x =us_start_x + us_slope_length
+    us_toe_start_y =top_start_y - us_slope_length / us_apron_skew
+    us_toe_end_y = top_start_y + bridge_width + us_slope_length / us_apron_skew
+
+    top_end_y = top_start_y + bridge_width
+    ds_slope_length = bank_hgt*ds_slope
+    ds_start_x = us_end_x + bank_width
+    ds_end_x = ds_start_x + ds_slope_length
+
+    ds_toe_start_y =top_start_y - ds_slope_length / ds_apron_skew
+    ds_toe_end_y = top_start_y + bridge_width + ds_slope_length / ds_apron_skew
+
+
+    N = len(x)
+    for i in range(N):
+
+       # Sloping Embankment Across Channel
+        if us_start_x < x[i] < us_end_x +0.1:   # For UPSLOPE on the Upstream FACE
+        #if 5.0 < x[i] < 10.1: # For a Range of X, and over a Range of Y based on X adjust Z
+            if us_toe_start_y +(x[i] - us_start_x)/us_apron_skew < y[i] < us_toe_end_y - (x[i] - us_start_x)/us_apron_skew:
+                #if  49.0+(x[i]-5.0)/5.0 <  y[i]  < 151.0 - (x[i]-5.0)/5.0: # Cut Out Base Segment for Culvert FACE
+                 z[i]=z[i] # Flat Apron
+                #z[i] += z[i] + (x[i] - us_start_x)/us_slope
+                #pass
+            else:
+               z[i] += z[i] + (x[i] - us_start_x)/us_slope    # Sloping Segment  U/S Face
+        if us_end_x < x[i] < ds_start_x + 0.1:
+           z[i] +=  z[i]+bank_hgt        # Flat Crest of Embankment
+        if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face
+            if  top_start_y-(x[i]-ds_start_x)/ds_apron_skew <  y[i]  < top_end_y + (x[i]-ds_start_x)/ds_apron_skew: # Cut Out Segment for Culvert FACE
+                 z[i]=z[i] # Flat Apron
+                #z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope
+                #pass
+            else:
+               z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope       # Sloping D/S Face
+           
+        
+
+    return z
+
+print 'Setting Quantities....'
+domain.set_quantity('elevation', topography)  # Use function for elevation
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+
+
+
+#------------------------------------------------------------------------------
+# Setup specialised forcing terms
+#------------------------------------------------------------------------------
+
+#------------------------------------------------------------------------------
+# Setup CULVERT INLETS and OUTLETS in Current Topography
+#------------------------------------------------------------------------------
+print 'Defining Structures'
+
+#  DEFINE CULVERT INLET AND OUTLETS
+
+
+#culvert0 = Culvert_operator(domain,
+#                            end_point0=[40.0, 75.0],
+#                            end_point1=[50.0, 75.0],
+#                            width=50.0,
+#                            depth=10.0,
+#                            apron=5.0,
+#                            verbose=False)
+
+
+#------------------------------------------------------------------------------
+# Setup culverts
+#------------------------------------------------------------------------------
+
+culverts = []
+number_of_culverts = 2
+for i in range(number_of_culverts):
+    culvert_width = 50.0/number_of_culverts
+    y = 100-i*culvert_width - culvert_width/2.0
+    ep0 = num.array([40.0, y])
+    ep1 = num.array([50.0, y])
+    
+    losses = {'inlet':0.5, 'outlet':1, 'bend':0, 'grate':0, 'pier': 0, 'other': 0}
+#    culverts.append(Boyd_pipe_operator(domain,
+#                            end_point0=ep0,
+#                            end_point1=ep1,
+#                            losses=losses,
+#                            diameter=1.5, #culvert_width, #3.658,
+#                            apron=6.0,
+#                            use_momentum_jet=True,
+#                            use_velocity_head=True,
+#                            manning=0.013,
+#                            logging=True,
+#                            label='culvert',
+#                            verbose=False))
+
+
+
+
+
+    Boyd_box_operator(domain,
+                      losses=losses,
+                      width=culvert_width,
+                      height=10.0,
+                      end_points=[ep0,ep1],
+                      apron=6.0,
+                      manning=0.013,
+                      label='box_culvert',
+                      logging=True,
+                      verbose=False)
+
+#losses = {'inlet':1, 'outlet':1, 'bend':1, 'grate':1, 'pier': 1, 'other': 1}
+#culvert2 = Culvert_operator(domain,
+                            #end_point0=[40.0, 62.5],
+                            #end_point1=[50.0, 62.5],
+                            #losses,
+                            #width=25.0,
+                            #depth=10.0,
+                            #apron=5.0,
+                            #manning=0.013,
+                            #verbose=False)
+
+
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+print 'Setting Boundary Conditions'
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
+
+Btus = anuga.Dirichlet_boundary([20.0, 0, 0])           # Outflow water at 10.0
+Btds = anuga.Dirichlet_boundary([19.0, 0, 0])           # Outflow water at 9.0
+domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
+
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 1, finaltime = 100):
+    print domain.timestepping_statistics()
+
+    domain.print_operator_timestepping_statistics()
+
Index: /trunk/anuga_core/demos/structures/run_update_special_condition.py
===================================================================
--- /trunk/anuga_core/demos/structures/run_update_special_condition.py	(revision 9507)
+++ /trunk/anuga_core/demos/structures/run_update_special_condition.py	(revision 9507)
@@ -0,0 +1,184 @@
+import os.path
+import sys
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.geometry.polygon_function import Polygon_function
+        
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+import anuga
+
+from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
+from anuga.structures.inlet_operator import Inlet_operator
+                            
+#from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
+     
+from math import pi, pow, sqrt
+
+import numpy as num
+
+
+
+
+
+"""test_that_culvert_runs_rating
+
+This test exercises the culvert and checks values outside rating curve
+are dealt with       
+"""
+
+path = get_pathname_from_package('anuga.culvert_flows')    
+
+length = 40.
+width = 15.
+
+dx = dy = 0.5          # Resolution: Length of subdivisions on both axes
+
+points, vertices, boundary = rectangular_cross(int(length/dx),
+                                               int(width/dy),
+                                               len1=length, 
+                                               len2=width)
+domain = anuga.Domain(points, vertices, boundary)   
+domain.set_name('Test_culvert')                 # Output name
+domain.set_default_order(2)
+#domain.set_beta(1.5)
+
+
+#----------------------------------------------------------------------
+# Setup initial conditions
+#----------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=-x/1000
+    
+    N = len(x)
+    for i in range(N):
+
+       # Sloping Embankment Across Channel
+        if 5.0 < x[i] < 10.1:
+            # Cut Out Segment for Culvert face                
+            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
+               z[i]=z[i]
+            else:
+               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+        if 10.0 < x[i] < 12.1:
+           z[i] +=  2.5                    # Flat Crest of Embankment
+        if 12.0 < x[i] < 14.5:
+            # Cut Out Segment for Culvert face                
+            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
+               z[i]=z[i]
+            else:
+               z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                   
+        
+    return z
+
+
+line = [[0.01, 5.0], [0.01, 10.0]]
+Q = 50.0
+
+import anuga
+
+inlet = anuga.structures.inlet.Inlet(domain, line)
+
+domain.special_inlet = inlet
+
+def my_update_special_conditions(domain):
+
+    #print 'inside my domain update_special condition'
+
+    domain.special_inlet.set_depths(5.0)
+
+import types
+domain.update_special_conditions = types.MethodType(my_update_special_conditions, domain)
+
+domain.set_quantity('elevation', topography) 
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+filename=os.path.join(path, 'example_rating_curve.csv')
+
+
+
+#Boyd_pipe_operator(domain,
+#                            end_point0=[9.0, 2.5],
+#                            end_point1=[13.0, 2.5],
+#                            losses=1.5,
+#                            diameter=1.5,
+#                            apron=5.0,
+#                            use_momentum_jet=True,
+#                            use_velocity_head=False,
+#                            manning=0.013,
+#                            verbose=False)
+
+
+line = [[0.0, 5.0], [0.0, 10.0]]
+Q = 50.0
+#Inlet_operator(domain, line, Q)
+
+
+
+
+
+
+
+##-----------------------------------------------------------------------
+## Setup boundary conditions
+##-----------------------------------------------------------------------
+
+## Inflow based on Flow Depth and Approaching Momentum
+Bi = anuga.Dirichlet_boundary([2.0, 0.0, 0.0])
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+#Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
+
+## Upstream and downstream conditions that will exceed the rating curve
+## I.e produce delta_h outside the range [0, 10] specified in the the 
+## file example_rating_curve.csv
+#Btus = anuga.Time_boundary(domain, \
+            #lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
+#Btds = anuga.Time_boundary(domain, \
+            #lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
+#domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+
+
+##-----------------------------------------------------------------------
+## Evolve system through time
+##-----------------------------------------------------------------------
+
+#min_delta_w = sys.maxint 
+#max_delta_w = -min_delta_w
+for t in domain.evolve(yieldstep = 1.0, finaltime = 30):
+    domain.write_time()
+
+    #if domain.get_time() > 150.5 and domain.get_time() < 151.5 :
+        #Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
+        #domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
+
+    #delta_w = culvert.inlet.stage - culvert.outlet.stage
+    
+    #if delta_w > max_delta_w: max_delta_w = delta_w
+    #if delta_w < min_delta_w: min_delta_w = delta_w
+
+    #print domain.volumetric_balance_statistics()
+
+
+
+
+
+## Check that extreme values in rating curve have been exceeded
+## so that we know that condition has been exercised
+#assert min_delta_w < 0
+#assert max_delta_w > 10        
+
+
+#os.remove('Test_culvert.sww')
Index: /trunk/anuga_core/demos/structures/run_wide_bridge.py
===================================================================
--- /trunk/anuga_core/demos/structures/run_wide_bridge.py	(revision 9507)
+++ /trunk/anuga_core/demos/structures/run_wide_bridge.py	(revision 9507)
@@ -0,0 +1,283 @@
+""" Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction
+
+This example includes a Model Topography that shows a TYPICAL Headwall Configuration
+
+The aim is to change the Culvert Routine to Model more precisely the abstraction
+from a vertical face.
+
+The inflow must include the impact of Approach velocity.
+Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style
+abstraction
+
+"""
+print 'Starting.... Importing Modules...'
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import anuga
+
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.shallow_water.shallow_water_domain import Domain
+from anuga.shallow_water.forcing import Rainfall, Inflow
+#from anuga.shallow_water.forcing import Reflective_boundary
+#from anuga.shallow_water.forcing import Dirichlet_boundary
+#from anuga.shallow_water.forcing import Transmissive_boundary, Time_boundary
+
+
+#from anuga.culvert_flows.culvert_routines import weir_orifice_channel_culvert_model
+from math import pi,pow,sqrt
+
+import numpy as num
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+print 'Setting up domain'
+
+length = 200. #x-Dir
+width = 200.  #y-dir
+
+dx = dy = 2.0          # Resolution: Length of subdivisions on both axes
+#dx = dy = .5           # Resolution: Length of subdivisions on both axes
+#dx = dy = .5           # Resolution: Length of subdivisions on both axes
+#dx = dy = .1           # 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)   
+domain.set_name('Test_WIDE_BRIDGE')                 # Output name
+#domain.set_default_order(2)
+#omain.H0 = 0.01
+#domain.tight_slope_limiters = 1
+
+domain.set_flow_algorithm('2_0')
+
+print 'Size', len(domain)
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+def topography(x, y):
+    """Set up a weir
+    
+    A culvert will connect either side
+    """
+    # General Slope of Topography
+    z=10.0-x/100.0  # % Longitudinal Slope
+    
+    #       NOW Add bits and Pieces to topography
+    bank_hgt=10.0
+    bridge_width = 50.0
+    bank_width = 10.0
+    
+    us_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
+    us_start_x = 10.0
+    top_start_y = 50.0
+    us_slope = 3.0  #Horiz : Vertic
+    ds_slope = 3.0
+    ds_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
+    centre_line_y= top_start_y+bridge_width/2.0
+
+    # CALCULATE PARAMETERS TO FORM THE EMBANKMENT
+    us_slope_length = bank_hgt*us_slope
+    us_end_x =us_start_x + us_slope_length
+    us_toe_start_y =top_start_y - us_slope_length / us_apron_skew
+    us_toe_end_y = top_start_y + bridge_width + us_slope_length / us_apron_skew
+
+    top_end_y = top_start_y + bridge_width
+    ds_slope_length = bank_hgt*ds_slope
+    ds_start_x = us_end_x + bank_width
+    ds_end_x = ds_start_x + ds_slope_length
+
+    ds_toe_start_y =top_start_y - ds_slope_length / ds_apron_skew
+    ds_toe_end_y = top_start_y + bridge_width + ds_slope_length / ds_apron_skew
+
+
+    N = len(x)
+    for i in range(N):
+
+       # Sloping Embankment Across Channel
+        if us_start_x < x[i] < us_end_x +0.1:   # For UPSLOPE on the Upstream FACE
+        #if 5.0 < x[i] < 10.1: # For a Range of X, and over a Range of Y based on X adjust Z
+            if us_toe_start_y +(x[i] - us_start_x)/us_apron_skew < y[i] < us_toe_end_y - (x[i] - us_start_x)/us_apron_skew:
+                #if  49.0+(x[i]-5.0)/5.0 <  y[i]  < 151.0 - (x[i]-5.0)/5.0: # Cut Out Base Segment for Culvert FACE
+                 z[i]=z[i] # Flat Apron
+                #z[i] += z[i] + (x[i] - us_start_x)/us_slope
+                #pass
+            else:
+               z[i] += z[i] + (x[i] - us_start_x)/us_slope    # Sloping Segment  U/S Face
+        if us_end_x < x[i] < ds_start_x + 0.1:
+           z[i] +=  z[i]+bank_hgt        # Flat Crest of Embankment
+        if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face
+            if  top_start_y-(x[i]-ds_start_x)/ds_apron_skew <  y[i]  < top_end_y + (x[i]-ds_start_x)/ds_apron_skew: # Cut Out Segment for Culvert FACE
+                 z[i]=z[i] # Flat Apron
+                #z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope
+                #pass
+            else:
+               z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope       # Sloping D/S Face
+           
+        
+
+    return z
+
+print 'Setting Quantities....'
+domain.set_quantity('elevation', topography)  # Use function for elevation
+domain.set_quantity('friction', 0.01)         # Constant friction 
+domain.set_quantity('stage',
+                    expression='elevation')   # Dry initial condition
+
+
+
+
+#------------------------------------------------------------------------------
+# Setup specialised forcing terms
+#------------------------------------------------------------------------------
+
+#------------------------------------------------------------------------------
+# Setup CULVERT INLETS and OUTLETS in Current Topography
+#------------------------------------------------------------------------------
+print 'DEFINING any Structures if Required'
+
+#  DEFINE CULVERT INLET AND OUTLETS
+
+"""
+culvert_rating = Culvert_flow(domain,
+    culvert_description_filename='example_rating_curve.csv',
+    end_point0=[0.0, 75.0], 
+    end_point1=[50.0, 75.0],
+    verbose=True)
+
+culvert_energy = Culvert_flow(domain,
+    label='Culvert No. 1',
+    description='This culvert is a test unit 4m diameter',   
+    end_point0=[40.0, 75.0], 
+    end_point1=[50.0, 75.0],
+    width=4.0,
+    culvert_routine=boyd_generalised_culvert_model,  #culvert_routine=weir_orifice_channel_culvert_model,     
+    number_of_barrels=1,
+    number_of_smoothing_steps=10,
+    #update_interval=0.25,
+    log_file=True,
+    discharge_hydrograph=True,
+    use_velocity_head=False,
+    use_momentum_jet=False,
+    verbose=True)
+
+domain.forcing_terms.append(culvert_energy)
+"""
+#from anuga.structures.boyd_box_operator import Boyd_box_operator
+#culvert0 = Culvert_operator(domain,
+#                            end_point0=[40.0, 75.0],
+#                            end_point1=[50.0, 75.0],
+#                            width=50.0,
+#                            depth=10.0,
+#                            apron=5.0,
+#                            verbose=False)
+
+
+#------------------------------------------------------------------------------
+# Setup culverts
+#------------------------------------------------------------------------------
+culverts = []
+number_of_culverts = 1 
+for i in range(number_of_culverts):
+    culvert_width = 50.0/number_of_culverts
+    y = 100-i*culvert_width - culvert_width/2.0
+    ep0 = num.array([37.0, y])
+    ep1 = num.array([53.0, y])
+    culverts.append(anuga.Boyd_box_operator(domain,
+                                            losses=1.5,
+                                            width=3.658,
+                                            height=3.658,
+                                            end_points=[ep0, ep1],
+                                            apron=0.5,
+                                            manning=0.013,
+                                            enquiry_gap=1.0,
+                                            description='bridge culvert',
+                                            verbose=False))
+
+#culvert2 = Culvert_operator(domain,
+#                            end_point0=[40.0, 62.5],
+#                            end_point1=[50.0, 62.5],
+#                            width=25.0,
+#                            depth=10.0,
+#                            apron=5.0,
+#                            manning=0.013,
+#                            verbose=False)
+
+
+
+"""
+culvert_energy = Culvert_flow(domain,
+    label='Culvert No. 1',
+    description='This culvert is a test unit 50m Wide by 5m High',   
+    end_point0=[40.0, 75.0],
+    end_point1=[50.0, 75.0],
+    width=50.0,depth=5.0,
+    culvert_routine=boyd_generalised_culvert_model,  #culvert_routine=weir_orifice_channel_culvert_model,     
+    number_of_barrels=1,
+    number_of_smoothing_steps=10,
+    #update_interval=0.25,
+    log_file=True,
+    discharge_hydrograph=True,
+    use_velocity_head=False,
+    use_momentum_jet=False,
+    verbose=True)
+
+domain.forcing_terms.append(culvert_energy)
+"""
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+print 'Setting Boundary Conditions'
+Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
+
+Bo = anuga.Dirichlet_boundary([-5.0, 0, 0])           # Outflow water at -5.0
+Bd = anuga.Dirichlet_boundary([0,0,0])                # Outflow water at 0.0
+
+#Btus = Time_boundary(domain, lambda t: [0.0+ 1.025*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0])
+#Btds = Time_boundary(domain, lambda t: [0.0+ 0.0075*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0])
+
+Btus = anuga.Dirichlet_boundary([20.0, 0, 0])           # Outflow water at 20
+Btds = anuga.Dirichlet_boundary([19.0, 0, 0])           # Outflow water at 19
+domain.set_boundary({'left': Btus, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 1, finaltime = 100):
+    print domain.timestepping_statistics()
+    print domain.volumetric_balance_statistics()
+    for i, culvert in enumerate(culverts):
+        print 'culvert: ', i
+        print culvert.timestepping_statistics()
+    
+
+    
+"""    
+#import sys; sys.exit() 
+# Profiling code
+import time
+t0 = time.time()
+   
+s = 'for t in domain.evolve(yieldstep = 0.1, finaltime = 300): domain.write_time()'
+
+import profile, pstats
+FN = 'profile.dat'
+
+profile.run(s, FN)
+    
+print 'That took %.2f seconds' %(time.time()-t0)
+
+S = pstats.Stats(FN)
+#S.sort_stats('time').print_stats(20)
+s = S.sort_stats('cumulative').print_stats(30)
+
+print s
+"""
