Index: /branches/numpy/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /branches/numpy/anuga/abstract_2d_finite_volumes/domain.py	(revision 6516)
+++ /branches/numpy/anuga/abstract_2d_finite_volumes/domain.py	(revision 6517)
@@ -319,4 +319,5 @@
         return self.mesh.get_boundary_polygon(*args, **kwargs)
 
+    # FIXME(Ole): This doesn't seem to be required
     def get_number_of_triangles_per_node(self, *args, **kwargs):
         return self.mesh.get_number_of_triangles_per_node(*args, **kwargs)
@@ -348,4 +349,7 @@
     def statistics(self, *args, **kwargs):
         return self.mesh.statistics(*args, **kwargs)
+        
+    def get_extent(self, *args, **kwargs):
+        return self.mesh.get_extent(*args, **kwargs)    
 
     ##
Index: /branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
===================================================================
--- /branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py	(revision 6516)
+++ /branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py	(revision 6517)
@@ -88,6 +88,13 @@
     # FIXME (Ole): We should rename f to function to be consistent with
     # Transmissive_Momentum_Set_Stage_Boundary (cf posting by rrraman)
-    def __init__(self, domain = None, f = None):
+    def __init__(self, domain=None,
+                 f=None,
+                 default_boundary=None,
+                 verbose=False):
         Boundary.__init__(self)
+        self.default_boundary = default_boundary
+        self.default_boundary_invoked = False    # Flag
+        self.domain = domain
+        self.verbose = verbose
 
         try:
@@ -122,5 +129,36 @@
     def evaluate(self, vol_id=None, edge_id=None):
         # FIXME (Ole): I think this should be get_time(), see ticket:306
-        return self.f(self.domain.time)
+        try:
+            res = self.f(self.domain.time)
+        except Modeltime_too_early, e:
+            raise Modeltime_too_early, e
+        except Modeltime_too_late, e:
+            if self.default_boundary is None:
+                raise Exception, e # Reraise exception
+            else:
+                # Pass control to default boundary
+                res = self.default_boundary.evaluate(vol_id, edge_id)
+                
+                # Ensure that result cannot be manipulated
+                # This is a real danger in case the 
+                # default_boundary is a Dirichlet type 
+                # for instance. 
+                res = res.copy() 
+                
+                if self.default_boundary_invoked is False:
+                    if self.verbose:                
+                        # Issue warning the first time
+                        msg = '%s' %str(e)
+                        msg += 'Instead I will use the default boundary: %s\n'\
+                            %str(self.default_boundary) 
+                        msg += 'Note: Further warnings will be supressed'
+                        print msg
+               
+                    # FIXME (Ole): Replace this crude flag with
+                    # Python's ability to print warnings only once.
+                    # See http://docs.python.org/lib/warning-filter.html
+                    self.default_boundary_invoked = True
+
+        return res
 
 
@@ -299,9 +337,10 @@
                     if self.default_boundary_invoked is False:
                         # Issue warning the first time
-                        msg = '%s' %str(e)
-                        msg += 'Instead I will use the default boundary: %s\n'\
-                            %str(self.default_boundary) 
-                        msg += 'Note: Further warnings will be supressed'
-                        warn(msg)
+                        if self.verbose:
+                            msg = '%s' %str(e)
+                            msg += 'Instead I will use the default boundary: %s\n'\
+                                %str(self.default_boundary) 
+                            msg += 'Note: Further warnings will be supressed'
+                            print msg
                    
                         # FIXME (Ole): Replace this crude flag with
Index: /branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py	(revision 6516)
+++ /branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py	(revision 6517)
@@ -886,5 +886,6 @@
         # a crash in fittng, so disabled it until I can investigate further
         # Sorry. 23 Jan 2009. Logged as ticket:314
-        if False:
+        #if True: # Test will fail (31 Jan 2009) 
+        if False: # Test will pass (31 Jan 2009)
             # Use mesh as defined by domain
             # This used to cause problems for caching due to quantities
Index: /branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py
===================================================================
--- /branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py	(revision 6516)
+++ /branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py	(revision 6517)
@@ -542,4 +542,5 @@
                         conserved_quantities =\
                         ['stage', 'xmomentum', 'ymomentum'])
+        domain.set_default_order(1)
         domain.check_integrity()
 
@@ -650,4 +651,5 @@
                         conserved_quantities =\
                         ['stage', 'xmomentum', 'ymomentum'])
+        domain.set_default_order(1)                        
         domain.check_integrity()
 
@@ -800,4 +802,48 @@
                              [ 11.0,  11.0,  11.0],
                              [ 11.0,  11.0,  11.0]])
+                             
+    def test_that_mesh_methods_exist(self):
+        """test_that_mesh_methods_exist
+        
+        Test that relavent mesh methods are made available in 
+        domain through composition
+        """
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+
+        # Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary)                             
+        
+        
+        domain.get_centroid_coordinates()
+        domain.get_radii()
+        domain.get_areas()
+        domain.get_area()
+        domain.get_vertex_coordinates()
+        domain.get_triangles()
+        domain.get_nodes()
+        domain.get_number_of_nodes()
+        domain.get_normal(0,0)
+        domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]])
+        domain.get_disconnected_triangles()
+        domain.get_boundary_tags()
+        domain.get_boundary_polygon()
+        #domain.get_number_of_triangles_per_node()
+        domain.get_triangles_and_vertices_per_node()
+        domain.get_interpolation_object()
+        domain.get_tagged_elements()
+        domain.get_lone_vertices()
+        domain.get_unique_vertices()
+        g = domain.get_georeference()
+        domain.set_georeference(g)
+        domain.build_tagged_elements_dictionary()
+        domain.statistics()
+        domain.get_extent()
+
+        
+        
 
 #-------------------------------------------------------------
Index: /branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py
===================================================================
--- /branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 6516)
+++ /branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 6517)
@@ -91,7 +91,13 @@
                   str(verts)))
         self.assert_(num.allclose(num.array([nodes_absolute[1],
-                                             nodes_absolute[0],
-                                             nodes_absolute[2]]),
-                                  verts), msg)
+                                     nodes_absolute[0],
+                                     nodes_absolute[2]]), verts))
+        verts = domain.get_vertex_coordinates(triangle_id=0,
+                                              absolute=True)       
+        self.assert_(num.allclose(num.array([nodes_absolute[1],
+                                     nodes_absolute[0],
+                                     nodes_absolute[2]]), verts))
+        
+        
 
     def test_get_vertex_coordinates_triangle_id(self):
Index: /branches/numpy/anuga/abstract_2d_finite_volumes/util.py
===================================================================
--- /branches/numpy/anuga/abstract_2d_finite_volumes/util.py	(revision 6516)
+++ /branches/numpy/anuga/abstract_2d_finite_volumes/util.py	(revision 6517)
@@ -925,4 +925,5 @@
                     verbose = False):   
         
+    # FIXME(Ole): Shouldn't print statements here be governed by verbose?
     assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
     
@@ -971,5 +972,6 @@
             raise msg
 
-        print 'swwfile', swwfile
+        if verbose:
+            print 'swwfile', swwfile
 
         # Extract parent dir name and use as label
@@ -2489,4 +2491,10 @@
                                  base_name=base,
                                  verbose=verbose)
+    #print 'sww files just after get_all_swwfiles()', sww_files
+    # fudge to get SWW files in 'correct' order, oldest on the left
+    sww_files.sort()
+
+    if verbose:
+        print 'sww files', sww_files
     
     #to make all the quantities lower case for file_function
@@ -2497,5 +2505,21 @@
 
     core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
-    
+    gauge_file = out_name
+
+    heading = [quantity for quantity in quantities]
+    heading.insert(0,'time')
+    heading.insert(1,'hours')
+
+    #create a list of csv writers for all the points and write header
+    points_writer = []
+    for point_i,point in enumerate(points):
+        points_writer.append(writer(file(dir_name + sep + gauge_file
+                                         + point_name[point_i] + '.csv', "wb")))
+        points_writer[point_i].writerow(heading)
+    
+    if verbose: print 'Writing csv files'
+
+    quake_offset_time = None
+
     for sww_file in sww_files:
         sww_file = join(dir_name, sww_file+'.sww')
Index: /branches/numpy/anuga/alpha_shape/alpha_shape.py
===================================================================
--- /branches/numpy/anuga/alpha_shape/alpha_shape.py	(revision 6516)
+++ /branches/numpy/anuga/alpha_shape/alpha_shape.py	(revision 6517)
@@ -289,10 +289,10 @@
                        (denom[k]< EPSILON and  denom[k] > -EPSILON)]
 
-        if num.any(denom == 0.0):
-            raise AlphaError
-
-        dx = num.divide(y31*dist21 - y21*dist31, denom)
-        dy = num.divide(x21*dist31 - x31*dist21, denom)
-
+        if num.alltrue(denom != 0.0):               
+            dx = num.divide(y31*dist21 - y21*dist31,denom)
+            dy = num.divide(x21*dist31 - x31*dist21,denom)
+        else:
+            raise  AlphaError
+            
         self.triradius = 0.5*num.sqrt(dx*dx + dy*dy)
         #print "triangle radii", self.triradius
Index: /branches/numpy/anuga/caching/test_caching.py
===================================================================
--- /branches/numpy/anuga/caching/test_caching.py	(revision 6516)
+++ /branches/numpy/anuga/caching/test_caching.py	(revision 6517)
@@ -386,9 +386,16 @@
 
         # Check that hash value of callable objects don't change
+        # FIXME (Ole): The hash values do appear to change when OS
+        # and/or dependencies are upgraded
         if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
           # 64 bit hash values
-          f1hash = 1914027059797211698
-          f2hash = 1914027059807087171
+          f1hash = 7079146893884768701
+          f2hash = -6995306676314913340
+
+          # Prior to cluster upgrades Feb 2009
+          #f1hash = 1914027059797211698 
+          #f2hash = 1914027059807087171
         else:
+          # 32 bit hash values
           f1hash = -758136387
           f2hash = -11221564     
Index: /branches/numpy/anuga/compile_all.py
===================================================================
--- /branches/numpy/anuga/compile_all.py	(revision 6516)
+++ /branches/numpy/anuga/compile_all.py	(revision 6517)
@@ -27,2 +27,4 @@
 #execfile('test_all.py')
     
+if sys.platform == 'win32':
+    raw_input('Press any key')
Index: /branches/numpy/anuga/config.py
===================================================================
--- /branches/numpy/anuga/config.py	(revision 6516)
+++ /branches/numpy/anuga/config.py	(revision 6517)
@@ -91,6 +91,5 @@
         
 # FIXME (Ole) Maybe get rid of order altogether and use beta_w
-# ... and isn't it about time we make the default 2?
-default_order = 1
+default_order = 2
 
 ################################################################################
Index: /branches/numpy/anuga/coordinate_transforms/redfearn.py
===================================================================
--- /branches/numpy/anuga/coordinate_transforms/redfearn.py	(revision 6516)
+++ /branches/numpy/anuga/coordinate_transforms/redfearn.py	(revision 6517)
@@ -37,5 +37,6 @@
     return sign*dd, mm, ss
 
-def redfearn(lat, lon, false_easting=None, false_northing=None, zone=None):
+def redfearn(lat, lon, false_easting=None, false_northing=None,
+             zone=None, central_meridian=None, scale_factor=None):
     """Compute UTM projection using Redfearn's formula
 
@@ -47,4 +48,8 @@
     If zone is specified reproject lat and long to specified zone instead of
     standard zone
+    If meridian is specified, reproject lat and lon to that instead of zone. In this case
+    zone will be set to -1 to indicate non-UTM projection
+
+    Note that zone and meridian cannot both be specifed
     """
 
@@ -57,5 +62,9 @@
     a = 6378137.0                       #Semi major axis
     inverse_flattening = 298.257222101  #1/f
-    K0 = 0.9996                         #Central scale factor    
+    if scale_factor is None:
+        K0 = 0.9996                         #Central scale factor    
+    else:
+        K0 = scale_factor
+    #print 'scale', K0
     zone_width = 6                      #Degrees
 
@@ -138,9 +147,17 @@
     m = term1 + term2 + term3 + term4 #OK
 
-    #Zone
+    if zone is not None and central_meridian is not None:
+        msg = 'You specified both zone and central_meridian. Provide only one of them'
+        raise Exception, msg
+    
+    # Zone
     if zone is None:
         zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
 
-    central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
+    # Central meridian
+    if central_meridian is None:
+        central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
+    else:
+        zone = -1
 
     omega = (lon-central_meridian)*pi/180 #Relative longitude (radians)
Index: /branches/numpy/anuga/coordinate_transforms/test_geo_reference.py
===================================================================
--- /branches/numpy/anuga/coordinate_transforms/test_geo_reference.py	(revision 6516)
+++ /branches/numpy/anuga/coordinate_transforms/test_geo_reference.py	(revision 6517)
@@ -235,24 +235,28 @@
         y = 3.0
         
-        g = Geo_reference(56, x, y)
-        points = [[3.0,34.0], [64.0,6.0]]
-        new_points = g.get_absolute(points)
-
-        self.failUnless(isinstance(new_points, list), 'failed')
-        self.failUnless(type(new_points) == type(points), 'failed')
-        for point, new_point in map(None, points, new_points):
-            self.failUnless(point[0]+x == new_point[0], 'failed')
-            self.failUnless(point[1]+y == new_point[1], 'failed')
-
-        # test with no supplied offsets
+        g = Geo_reference(56,x,y)
+        lofl = [[3.0,34.0], [64.0,6.0]]
+        new_lofl = g.get_absolute(lofl)
+        #print "lofl",lofl 
+        #print "new_lofl",new_lofl
+
+        self.failUnless(type(new_lofl) == types.ListType, ' failed')
+        self.failUnless(type(new_lofl) == type(lofl), ' failed')
+        for point,new_point in map(None,lofl,new_lofl):
+            self.failUnless(point[0]+x==new_point[0], ' failed')
+            self.failUnless(point[1]+y==new_point[1], ' failed')
+
+            
         g = Geo_reference()
-        points = [[3.0,34.0], [64.0,6.0]]
-        new_points = g.get_absolute(points)
-
-        self.failUnless(isinstance(new_points, list), 'failed')
-        self.failUnless(type(new_points) == type(points), 'failed')
-        for point, new_point in map(None, points, new_points):
-            self.failUnless(point[0] == new_point[0], 'failed')
-            self.failUnless(point[1] == new_point[1], 'failed')
+        lofl = [[3.0,34.0], [64.0,6.0]]
+        new_lofl = g.get_absolute(lofl)
+        #print "lofl",lofl 
+        #print "new_lofl",new_lofl
+
+        self.failUnless(type(new_lofl) == types.ListType, ' failed')
+        self.failUnless(type(new_lofl) == type(lofl), ' failed')
+        for point,new_point in map(None,lofl,new_lofl):
+            self.failUnless(point[0]==new_point[0], ' failed')
+            self.failUnless(point[1]==new_point[1], ' failed')
             
         # test that calling get_absolute twice does the right thing
Index: /branches/numpy/anuga/coordinate_transforms/test_redfearn.py
===================================================================
--- /branches/numpy/anuga/coordinate_transforms/test_redfearn.py	(revision 6516)
+++ /branches/numpy/anuga/coordinate_transforms/test_redfearn.py	(revision 6517)
@@ -122,7 +122,14 @@
 
     def test_UTM_6_nonstandard_projection(self):
-        #Test 6 (Geraldton, WA)
-
-        #First test native projection (zone 50)
+        """test_UTM_6_nonstandard_projection
+
+        Test that projections can be forced to
+        use other than native zone.
+
+        Data is from Geraldton, WA
+        """
+        
+
+        # First test native projection (zone 50)
         zone, easting, northing = redfearn(-29.233299999,114.05)
 
@@ -267,4 +274,108 @@
     #    #assert allclose(northing, 6181725.1724276)
 
+    
+    def Xtest_nonstandard_meridian_coinciding_with_native(self):
+        """test_nonstandard_meridian_coinciding_with_native
+
+        This test will verify that redfearn can be used to project
+        points using an arbitrary central meridian that happens to 
+        coincide with the standard meridian at the center of a UTM zone.
+        This is a preliminary test before testing this functionality
+        with a truly arbitrary non-standard meridian.
+        """
+
+        # The file projection_test_points_z53.csv contains 10 points
+        # which straddle the boundary between UTM zones 53 and 54.
+        # They have been projected to zone 53 irrespective of where they
+        # belong.
+
+        path = get_pathname_from_package('anuga.coordinate_transforms')
+        
+        for forced_zone in [53, 54]:
+        
+            datafile = join(path, 'projection_test_points_z%d.csv' % forced_zone)
+            fid = open(datafile)
+
+            for line in fid.readlines()[1:]:
+                fields = line.strip().split(',')
+                
+                lon = float(fields[1])
+                lat = float(fields[2])
+                x = float(fields[3])
+                y = float(fields[4])            
+
+                zone, easting, northing = redfearn(lat, lon,
+                                                   zone=forced_zone)
+                
+                print
+                print 'Lat', lat
+                print 'Lon', lon
+                print 'Zone', zone
+                print 'Ref x', x, 'Computed x', easting, 'Close enough:', num.allclose(x, easting)
+                print 'Ref y', y, 'Computed y', northing, 'Close enough:', num.allclose(y, northing)
+                
+                # Check calculation
+                assert zone == forced_zone
+                print 
+                #assert num.allclose(x, easting)
+                #assert num.allclose(y, northing)
+
+    
+    
+    
+    def Xtest_nonstandard_meridian(self):
+        """test_nonstandard_meridian
+
+        This test will verify that redfearn can be used to project
+        points using an arbitrary central meridian.
+        """
+
+        # The file projection_test_points.csv contains 10 points
+        # which straddle the boundary between UTM zones 53 and 54.
+        # They have been projected using a central meridian of 137.5
+        # degrees (the boundary is 138 so it is pretty much right
+        # in the middle of zones 53 and 54).
+
+        path = get_pathname_from_package('anuga.coordinate_transforms')
+        datafile = join(path, 'projection_test_points.csv')
+        fid = open(datafile)
+
+        for line in fid.readlines()[1:]:
+            fields = line.strip().split(',')
+
+            lon = float(fields[1])
+            lat = float(fields[2])
+            x = float(fields[3])
+            y = float(fields[4])            
+
+            zone, easting, northing = redfearn(lat, lon,
+                                               central_meridian=137.5,
+                                               scale_factor=0.9996)
+
+            print
+            print 'Lat', lat
+            print 'Lon', lon
+            print 'Zone', zone
+            print 'Ref x', x, 'Computed x', easting, 'Close enough:', num.allclose(x, easting)
+            print 'Ref y', y, 'Computed y', northing, 'Close enough:', num.allclose(y, northing)
+
+            # Check calculation
+            assert zone == -1 # Indicates non UTM projection
+            print 
+            #assert num.allclose(x, easting)
+            #assert num.allclose(y, northing)
+
+        # Test that zone and meridian can't both be specified
+        try:
+            zone, easting, northing = redfearn(lat, lon,
+                                               zone=50, 
+                                               central_meridian=137.5)
+        except:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise Exception, msg
+
+            
     def test_convert_lats_longs(self):
 
Index: /branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py
===================================================================
--- /branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py	(revision 6516)
+++ /branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py	(revision 6517)
@@ -123,4 +123,6 @@
                        number_of_barrels=1,
                        update_interval=2,
+                       log_file=True,
+                       discharge_hydrograph=True,
                        verbose=True)
 
Index: /branches/numpy/anuga/culvert_flows/culvert_class.py
===================================================================
--- /branches/numpy/anuga/culvert_flows/culvert_class.py	(revision 6516)
+++ /branches/numpy/anuga/culvert_flows/culvert_class.py	(revision 6517)
@@ -204,4 +204,6 @@
             log_to_file(self.log_filename, description)
             log_to_file(self.log_filename, self.culvert_type)        
+        else:
+            self.log_filename = None
 
 
@@ -296,5 +298,5 @@
 
         # Print some diagnostics to log if requested
-        if hasattr(self, 'log_filename'):
+        if self.log_filename is not None:
             s = 'Culvert Effective Length = %.2f m' %(self.length)
             log_to_file(self.log_filename, s)
@@ -324,5 +326,6 @@
                 update = True
 
-        if hasattr(self, 'log_filename'):            
+            
+        if self.log_filename is not None:        
             s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
             log_to_file(self.log_filename, s)
@@ -439,5 +442,6 @@
             if self.verbose is True:
                 print msg
-            if hasattr(self, 'log_filename'):                    
+                
+            if self.log_filename is not None:                
                 log_to_file(self.log_filename, msg)
         
@@ -524,5 +528,5 @@
                 print '%.2fs - WARNING: Flow is running uphill.' % time
             
-        if hasattr(self, 'log_filename'):
+        if self.log_filename is not None:
             s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
                 %(time, self.inlet.stage, self.outlet.stage)
@@ -573,10 +577,21 @@
                     msg += 'for culvert "%s"' % self.label 
                     
-                    if hasattr(self, 'log_filename'):                    
+                    if self.log_filename is not None:                    
                         log_to_file(self.log_filename, msg)
             else:
                 # User culvert routine
                 Q, barrel_velocity, culvert_outlet_depth =\
-                    self.culvert_routine(self, delta_total_energy, g)
+                    self.culvert_routine(inlet.depth,
+                                         outlet.depth,
+                                         inlet.specific_energy, 
+                                         delta_total_energy, 
+                                         g,
+                                         culvert_length=self.length,
+                                         culvert_width=self.width,
+                                         culvert_height=self.height,
+                                         culvert_type=self.culvert_type,
+                                         manning=self.manning,
+                                         sum_loss=self.sum_loss,
+                                         log_filename=self.log_filename)
             
             
@@ -601,5 +616,5 @@
 
                             
-    
+# OBSOLETE (Except for momentum jet in Culvert_flow_energy)    
 class Culvert_flow_rating:
     """Culvert flow - transfer water from one hole to another
Index: /branches/numpy/anuga/culvert_flows/culvert_routines.py
===================================================================
--- /branches/numpy/anuga/culvert_flows/culvert_routines.py	(revision 6516)
+++ /branches/numpy/anuga/culvert_flows/culvert_routines.py	(revision 6517)
@@ -8,31 +8,8 @@
 
 #NOTE:
-# Inlet control:  Delta_Et > Es at the inlet
-# Outlet control: Delta_Et < Es at the inlet
-# where Et is total energy (w + 0.5*v^2/g) and
-# Es is the specific energy (h + 0.5*v^2/g)
-
-
-
-#   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
-# 
-# The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
-# Flow Is Removed at a Rate of INFLOW 
-# Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
-#
-# This SHould be changed to a Vertical Opening Both BOX and Circular
-# There will be several Culvert Routines such as:
-# CULVERT_Boyd_Channel
-# CULVERT_Orifice_and_Weir
-# CULVERT_Simple_FLOOR
-# CULVERT_Simple_WALL
-# CULVERT_Eqn_Floor
-# CULVERT_Eqn_Wall
-# CULVERT_Tab_Floor
-# CULVERT_Tab_Wall
-# BRIDGES.....
-# NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
-
-#  COULD USE EPA SWMM Model !!!!
+# Inlet control:  Delta_total_energy > inlet_specific_energy
+# Outlet control: Delta_total_energy < inlet_specific_energy
+# where total energy is (w + 0.5*v^2/g) and
+# specific energy is (h + 0.5*v^2/g)
 
 
@@ -40,11 +17,27 @@
 
 
-def boyd_generalised_culvert_model(culvert, delta_total_energy, g):
-
-    """Boyd's generalisation of the US department of transportation culvert model
-        # == The quantity of flow passing through a culvert is controlled by many factors
-        # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
-        # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
-        # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled 
+def boyd_generalised_culvert_model(inlet_depth, 
+                                   outlet_depth,
+                                   inlet_specific_energy, 
+                                   delta_total_energy, 
+                                   g,
+                                   culvert_length=0.0,
+                                   culvert_width=0.0,
+                                   culvert_height=0.0,
+                                   culvert_type='box',
+                                   manning=0.0,
+                                   sum_loss=0.0,
+                                   log_filename=None):
+
+    """Boyd's generalisation of the US department of transportation culvert 
+    model
+      
+    The quantity of flow passing through a culvert is controlled by many factors
+    It could be that the culvert is controlled by the inlet only, with it 
+    being unsubmerged this is effectively equivalent to the weir Equation
+    Else the culvert could be controlled by the inlet, with it being 
+    submerged, this is effectively the Orifice Equation
+    Else it may be controlled by down stream conditions where depending on 
+    the down stream depth, the momentum in the culvert etc. flow is controlled 
     """
 
@@ -53,55 +46,33 @@
     from anuga.utilities.numerical_tools import safe_acos as acos
 
-    inlet = culvert.inlet
-    outlet = culvert.outlet        
-    Q_outlet_tailwater = 0.0
-    inlet.rate = 0.0
-    outlet.rate = 0.0
-    Q_inlet_unsubmerged = 0.0
-    Q_inlet_submerged = 0.0
-    Q_outlet_critical_depth = 0.0
-
-    if hasattr(culvert, 'log_filename'):                                
-        log_filename = culvert.log_filename
-
-    manning = culvert.manning
-    sum_loss = culvert.sum_loss
-    length = culvert.length
-
-    if inlet.depth > 0.01:
-
-        E = inlet.specific_energy
-
-        if hasattr(culvert, 'log_filename'):                            
-            s = 'Specific energy  = %f m' %E
+    if inlet_depth > 0.01:
+        # Water has risen above inlet
+        
+        if log_filename is not None:                            
+            s = 'Specific energy  = %f m' % inlet_specific_energy
             log_to_file(log_filename, s)
         
-        msg = 'Specific energy is negative'
-        assert E >= 0.0, msg
+        msg = 'Specific energy at inlet is negative'
+        assert inlet_specific_energy >= 0.0, msg
                       
-        
-        # Water has risen above inlet
-        if culvert.culvert_type == 'circle':
-            # Round culvert
+        if culvert_type == 'circle':
+            # Round culvert (use width as diameter)
+            diameter = culvert_width            
 
             # Calculate flows for inlet control            
-            diameter = culvert.diameter
-
-            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*E**1.63    # Inlet Ctrl Inlet Unsubmerged 
-            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**0.63      # Inlet Ctrl Inlet Submerged
-
-            if hasattr(culvert, 'log_filename'):
-                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
+            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged 
+            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
+
+            if log_filename is not None:                                        
+                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
                 log_to_file(log_filename, s)
 
-            case = ''
+            # FIXME(Ole): Are these functions really for inlet control?                    
             if Q_inlet_unsubmerged < Q_inlet_submerged:
                 Q = Q_inlet_unsubmerged
 
-                alpha = acos(1 - inlet.depth/diameter)
+                alpha = acos(1 - inlet_depth/diameter)
                 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
-                outlet_culvert_depth = inlet.depth
-                width = diameter*sin(alpha)
-                #perimeter = alpha*diameter                
+                outlet_culvert_depth = inlet_depth
                 case = 'Inlet unsubmerged'
             else:    
@@ -109,47 +80,42 @@
                 flow_area = (diameter/2)**2 * pi
                 outlet_culvert_depth = diameter
-                width = diameter
-                #perimeter = diameter
                 case = 'Inlet submerged'                    
 
 
 
-            if delta_total_energy < E:
+            if delta_total_energy < inlet_specific_energy:
                 # Calculate flows for outlet control
                 # Determine the depth at the outlet relative to the depth of flow in the Culvert
 
-                if outlet.depth > diameter:        # The Outlet is Submerged
+                if outlet_depth > diameter:        # The Outlet is Submerged
                     outlet_culvert_depth=diameter
                     flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
                     perimeter = diameter * pi
-                    width = diameter
                     case = 'Outlet submerged'
-                elif outlet.depth==0.0: 
-                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
-                    alpha = acos(1 - inlet.depth/diameter)
+                elif outlet_depth==0.0: 
+                    outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth
+                    alpha = acos(1 - inlet_depth/diameter)
                     flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
                     perimeter = alpha*diameter
-                    width = diameter*sin(alpha)                    
 
                     case = 'Outlet depth is zero'                        
                 else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
-                    outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
-                    alpha = acos(1 - outlet.depth/diameter)
+                    outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
+                    alpha = acos(1 - outlet_depth/diameter)
                     flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
                     perimeter = alpha*diameter
-                    width = diameter*sin(alpha)                    
                     case = 'Outlet is open channel flow'
 
                 hyd_rad = flow_area/perimeter
                 
-                if hasattr(culvert, 'log_filename'):
+                if log_filename is not None:
                     s = 'hydraulic radius at outlet = %f' %hyd_rad
                     log_to_file(log_filename, s)
 
                 # Outlet control velocity using tail water
-                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
+                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
                 Q_outlet_tailwater = flow_area * culvert_velocity
                 
-                if hasattr(culvert, 'log_filename'):
+                if log_filename is not None:                
                     s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
                     log_to_file(log_filename, s)
@@ -165,20 +131,19 @@
 
             # Calculate flows for inlet control
-            height = culvert.height
-            width = culvert.width            
+            height = culvert_height
+            width = culvert_width            
             
-            Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
-            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
-
-            if hasattr(culvert, 'log_filename'):
+            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
+            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
+
+            if log_filename is not None:                            
                 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
                 log_to_file(log_filename, s)
 
-            case = ''
+            # FIXME(Ole): Are these functions really for inlet control?    
             if Q_inlet_unsubmerged < Q_inlet_submerged:
                 Q = Q_inlet_unsubmerged
-                flow_area = width*inlet.depth
-                outlet_culvert_depth = inlet.depth
-                #perimeter=(width+2.0*inlet.depth)                
+                flow_area = width*inlet_depth
+                outlet_culvert_depth = inlet_depth
                 case = 'Inlet unsubmerged'
             else:    
@@ -189,36 +154,36 @@
                 case = 'Inlet submerged'                    
 
-            if delta_total_energy < E:
+            if delta_total_energy < inlet_specific_energy:
                 # Calculate flows for outlet control
                 # Determine the depth at the outlet relative to the depth of flow in the Culvert
 
-                if outlet.depth > height:        # The Outlet is Submerged
+                if outlet_depth > height:        # The Outlet is Submerged
                     outlet_culvert_depth=height
                     flow_area=width*height       # Cross sectional area of flow in the culvert
                     perimeter=2.0*(width+height)
                     case = 'Outlet submerged'
-                elif outlet.depth==0.0: 
-                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
-                    flow_area=width*inlet.depth
-                    perimeter=(width+2.0*inlet.depth)
+                elif outlet_depth==0.0: 
+                    outlet_culvert_depth=inlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
+                    flow_area=width*inlet_depth
+                    perimeter=(width+2.0*inlet_depth)
                     case = 'Outlet depth is zero'                        
                 else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
-                    outlet_culvert_depth=outlet.depth
-                    flow_area=width*outlet.depth
-                    perimeter=(width+2.0*outlet.depth)
+                    outlet_culvert_depth=outlet_depth
+                    flow_area=width*outlet_depth
+                    perimeter=(width+2.0*outlet_depth)
                     case = 'Outlet is open channel flow'
 
                 hyd_rad = flow_area/perimeter
                 
-                if hasattr(culvert, 'log_filename'):                
-                    s = 'hydraulic radius at outlet = %f' %hyd_rad
+                if log_filename is not None:                                
+                    s = 'hydraulic radius at outlet = %f' % hyd_rad
                     log_to_file(log_filename, s)
 
                 # Outlet control velocity using tail water
-                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
+                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
                 Q_outlet_tailwater = flow_area * culvert_velocity
 
-                if hasattr(culvert, 'log_filename'):                            
-                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
+                if log_filename is not None:                            
+                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
                     log_to_file(log_filename, s)
                 Q = min(Q, Q_outlet_tailwater)
@@ -228,20 +193,15 @@
 
         # Common code for circle and square geometries
-        if hasattr(culvert, 'log_filename'):                                    
-            log_to_file(log_filename, 'Case: "%s"' %case)
+        if log_filename is not None:
+            log_to_file(log_filename, 'Case: "%s"' % case)
             
-        flow_rate_control=Q
-
-        if hasattr(culvert, 'log_filename'):                                    
-            s = 'Flow Rate Control = %f' %flow_rate_control
+        if log_filename is not None:                        
+            s = 'Flow Rate Control = %f' % Q
             log_to_file(log_filename, s)
 
-        inlet.rate = -flow_rate_control
-        outlet.rate = flow_rate_control                
-            
-        culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
-        
-        if hasattr(culvert, 'log_filename'):                            
-            s = 'Froude in Culvert = %f' %culv_froude
+        culv_froude=sqrt(Q**2*width/(g*flow_area**3))
+        
+        if log_filename is not None:                            
+            s = 'Froude in Culvert = %f' % culv_froude
             log_to_file(log_filename, s)
 
@@ -250,5 +210,5 @@
 
 
-    else: #inlet.depth < 0.01:
+    else: # inlet_depth < 0.01:
         Q = barrel_velocity = outlet_culvert_depth = 0.0
         
Index: /branches/numpy/anuga/fit_interpolate/fit.py
===================================================================
--- /branches/numpy/anuga/fit_interpolate/fit.py	(revision 6516)
+++ /branches/numpy/anuga/fit_interpolate/fit.py	(revision 6517)
@@ -310,4 +310,5 @@
                 assert flag is True, msg                
                 
+                # FIXME(Ole): This is the message referred to in ticket:314
                 minx = min(self.mesh_boundary_polygon[:,0])
                 maxx = max(self.mesh_boundary_polygon[:,0])                
@@ -317,4 +318,5 @@
                 msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
                     % (minx, maxx, miny, maxy)
+                raise RuntimeError, msg
                 
 
@@ -375,4 +377,12 @@
                 self.build_fit_subset(points, z, verbose=verbose)
 
+                # FIXME(Ole): I thought this test would make sense here
+                # See test_fitting_example_that_crashed_2 in test_shallow_water_domain.py
+                # Committed 11 March 2009
+                msg = 'Matrix AtA was not built'
+                assert self.AtA is not None, msg
+                
+                #print 'Matrix was built OK'
+
                 
             point_coordinates = None
@@ -382,8 +392,7 @@
         if point_coordinates is None:
             if verbose: print 'Warning: no data points in fit'
-            #assert self.AtA != None, 'no interpolation matrix'
-            #assert self.Atz != None
-            assert not self.AtA is None, 'no interpolation matrix'
-            assert not self.Atz is None
+            msg = 'No interpolation matrix.'
+            assert self.AtA is not None, msg
+            assert self.Atz is not None
             
             # FIXME (DSG) - do  a message
Index: /branches/numpy/anuga/geospatial_data/geospatial_data.py
===================================================================
--- /branches/numpy/anuga/geospatial_data/geospatial_data.py	(revision 6516)
+++ /branches/numpy/anuga/geospatial_data/geospatial_data.py	(revision 6517)
@@ -1806,4 +1806,17 @@
     points = G_small.get_data_points()
 
+    #FIXME: Remove points outside boundary polygon
+#    print 'new point',len(points)
+#
+#    new_points=[]
+#    new_points=array([],typecode=Float)
+#    new_points=resize(new_points,(len(points),2))
+#    print "BOUNDARY", boundary_poly
+#    for i,point in enumerate(points):
+#        if is_inside_polygon(point,boundary_poly, verbose=True):
+#            new_points[i] = point
+#            print"WOW",i,new_points[i]
+#    points = new_points
+
     if verbose: print "Number of points in sample to compare: ", len(points)
 
Index: /branches/numpy/anuga/shallow_water/__init__.py
===================================================================
--- /branches/numpy/anuga/shallow_water/__init__.py	(revision 6516)
+++ /branches/numpy/anuga/shallow_water/__init__.py	(revision 6517)
@@ -12,5 +12,6 @@
      Transmissive_momentum_set_stage_boundary,\
      Dirichlet_discharge_boundary,\
-     Field_boundary
+     Field_boundary,\
+     Transmissive_stage_zero_momentum_boundary
 
 
Index: /branches/numpy/anuga/shallow_water/data_manager.py
===================================================================
--- /branches/numpy/anuga/shallow_water/data_manager.py	(revision 6516)
+++ /branches/numpy/anuga/shallow_water/data_manager.py	(revision 6517)
@@ -5863,5 +5863,6 @@
         # This is being used to seperate one number from a list.
         # what it is actually doing is sorting lists from numeric arrays.
-        if type(times) is list or isinstance(times, num.ndarray):
+        #if type(times) is list or isinstance(times, num.ndarray):
+        if isinstance(times, (list, num.ndarray)):
             number_of_times = len(times)
             times = ensure_numeric(times)
@@ -5930,6 +5931,7 @@
             #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
 
-        if type(times) is list or isinstance(times, num.ndarray):
-            outfile.variables['time'][:] = times    #Store time relative
+        #if type(times) is list or isinstance(times, num.ndarray):
+        if isinstance(times, (list, num.ndarray)):
+            outfile.variables['time'][:] = times    # Store time relative
 
         if verbose:
@@ -6265,5 +6267,6 @@
         # This is being used to seperate one number from a list.
         # what it is actually doing is sorting lists from numeric arrays.
-        if type(times) is list or isinstance(times, num.ndarray):
+        #if type(times) is list or isinstance(times, num.ndarray):
+        if isinstance(times, (list, num.ndarray)):
             number_of_times = len(times)
             times = ensure_numeric(times)
@@ -6312,5 +6315,6 @@
             outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
 
-        if type(times) is list or isinstance(times, num.ndarray):
+        #if type(times) is list or isinstance(times, num.ndarray):
+        if isinstance(times, (list, num.ndarray)):
             outfile.variables['time'][:] = times    #Store time relative
 
Index: /branches/numpy/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- /branches/numpy/anuga/shallow_water/test_shallow_water_domain.py	(revision 6516)
+++ /branches/numpy/anuga/shallow_water/test_shallow_water_domain.py	(revision 6517)
@@ -1672,6 +1672,7 @@
         # Setup computational domain
         #-----------------------------------------------------------------
-        points, vertices, boundary = rectangular_cross(10, 10)    # Basic mesh
-        domain = Domain(points, vertices, boundary)    # Create domain
+        points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
+        domain = Domain(points, vertices, boundary) # Create domain
+        domain.set_default_order(1)        
         domain.set_quantities_to_be_stored(None)
         domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
Index: /branches/numpy/anuga/test_all.py
===================================================================
--- /branches/numpy/anuga/test_all.py	(revision 6516)
+++ /branches/numpy/anuga/test_all.py	(revision 6517)
@@ -215,2 +215,5 @@
         #os.remove(filename)
 
+    
+    if sys.platform == 'win32':
+        raw_input('Press any key')
