Index: /trunk/anuga_core/source/anuga/structures/boyd_box_culvert.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/boyd_box_culvert.py	(revision 7983)
+++ /trunk/anuga_core/source/anuga/structures/boyd_box_culvert.py	(revision 7984)
@@ -15,9 +15,11 @@
                  domain,
                  end_points,
-                 width=None,
-                 height=None,
+                 width,
+                 height,
+                 apron,
+                 enquiry_gap,
                  verbose=False):
 
-        culvert.Culvert.__init__(self, domain, end_points, width, height, verbose)
+        culvert.Culvert.__init__(self, domain, end_points, width, height, apron, enquiry_gap, verbose)
 
         self.routine = boyd_box_routine.Boyd_box_routine(self)
Index: /trunk/anuga_core/source/anuga/structures/boyd_box_culvert_routine.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/boyd_box_culvert_routine.py	(revision 7983)
+++ /trunk/anuga_core/source/anuga/structures/boyd_box_culvert_routine.py	(revision 7984)
@@ -70,8 +70,8 @@
         local_debug ='false'
         
-        if self.inflow.get_average_height() > 0.01: #this value was 0.01:
+        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01:
             if local_debug =='true':
                 log.critical('Specific E & Deltat Tot E = %s, %s'
-                             % (str(self.inflow.get_average_specific_energy()),
+                             % (str(self.inflow.get_enquiry_specific_energy()),
                                 str(self.delta_total_energy)))
                 log.critical('culvert type = %s' % str(culvert_type))
@@ -79,9 +79,9 @@
 
             if self.log_filename is not None:
-                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
+                s = 'Specific energy  = %f m' % self.inflow.get_enquiry_specific_energy()
                 log_to_file(self.log_filename, s)
 
             msg = 'Specific energy at inlet is negative'
-            assert self.inflow.get_average_specific_energy() >= 0.0, msg
+            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
 
             height = self.culvert_height
@@ -89,6 +89,6 @@
             flow_width = self.culvert_width
 
-            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
-            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_average_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
+            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_enquiry_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
+            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_enquiry_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
 
             # FIXME(Ole): Are these functions really for inlet control?
@@ -120,9 +120,9 @@
                 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
 
-            if self.delta_total_energy < self.inflow.get_average_specific_energy():
+            if self.delta_total_energy < self.inflow.get_enquiry_specific_energy():
                 # Calculate flows for outlet control
 
                 # Determine the depth at the outlet relative to the depth of flow in the Culvert
-                if self.outflow.get_average_height() > height:        # The Outlet is Submerged
+                if self.outflow.get_enquiry_height() > height:        # The Outlet is Submerged
                     outlet_culvert_depth=height
                     flow_area=width*height       # Cross sectional area of flow in the culvert
@@ -172,5 +172,5 @@
         # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
 
-        else: # self.inflow.get_average_height() < 0.01:
+        else: # self.inflow.get_enquiry_height() < 0.01:
             Q = barrel_velocity = outlet_culvert_depth = 0.0
 
Index: /trunk/anuga_core/source/anuga/structures/boyd_box_routine.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/boyd_box_routine.py	(revision 7983)
+++ /trunk/anuga_core/source/anuga/structures/boyd_box_routine.py	(revision 7984)
@@ -70,8 +70,8 @@
         local_debug ='false'
         
-        if self.inflow.get_average_height() > 0.01: #this value was 0.01:
+        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01:
             if local_debug =='true':
                 log.critical('Specific E & Deltat Tot E = %s, %s'
-                             % (str(self.inflow.get_average_specific_energy()),
+                             % (str(self.inflow.get_enquiry_specific_energy()),
                                 str(self.delta_total_energy)))
                 log.critical('culvert type = %s' % str(culvert_type))
@@ -79,9 +79,9 @@
 
             if self.log_filename is not None:
-                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
+                s = 'Specific energy  = %f m' % self.inflow.get_enquiry_specific_energy()
                 log_to_file(self.log_filename, s)
 
             msg = 'Specific energy at inlet is negative'
-            assert self.inflow.get_average_specific_energy() >= 0.0, msg
+            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
 
             height = self.culvert_height
@@ -89,6 +89,6 @@
             flow_width = self.culvert_width
 
-            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
-            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_average_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
+            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_enquiry_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
+            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_enquiry_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
 
             # FIXME(Ole): Are these functions really for inlet control?
@@ -120,9 +120,9 @@
                 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
 
-            if self.delta_total_energy < self.inflow.get_average_specific_energy():
+            if self.delta_total_energy < self.inflow.get_enquiry_specific_energy():
                 # Calculate flows for outlet control
 
                 # Determine the depth at the outlet relative to the depth of flow in the Culvert
-                if self.outflow.get_average_height() > height:        # The Outlet is Submerged
+                if self.outflow.get_enquiry_height() > height:        # The Outlet is Submerged
                     outlet_culvert_depth=height
                     flow_area=width*height       # Cross sectional area of flow in the culvert
@@ -172,5 +172,5 @@
         # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
 
-        else: # self.inflow.get_average_height() < 0.01:
+        else: # self.inflow.get_enquiry_height() < 0.01:
             Q = barrel_velocity = outlet_culvert_depth = 0.0
 
Index: /trunk/anuga_core/source/anuga/structures/culvert.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/culvert.py	(revision 7983)
+++ /trunk/anuga_core/source/anuga/structures/culvert.py	(revision 7984)
@@ -20,17 +20,18 @@
                  domain,
                  end_points, 
-                 width=None,
-                 height=None,
-                 verbose=False):
+                 width,
+                 height,
+                 apron,
+                 enquiry_gap,
+                 verbose):
         
         # Input check
         
         self.domain = domain
-
         self.end_points = end_points
- 
-        self.width = width
+        self.width  = width
         self.height = height
-        
+        self.apron  = apron
+        self.enquiry_gap = enquiry_gap
         self.verbose=verbose
        
@@ -42,10 +43,12 @@
         self.inlets = []
         polygon0 = self.inlet_polygons[0]
+        ep0 = self.inlet_equiry_pts[0]
         outward_vector0 = self.culvert_vector
-        self.inlets.append(inlet.Inlet(self.domain, polygon0, outward_vector0))
+        self.inlets.append(inlet.Inlet(self.domain, polygon0, ep0, outward_vector0))
 
         polygon1 = self.inlet_polygons[1]
+        ep1 = self.inlet_equiry_pts[1]
         outward_vector1  = - self.culvert_vector
-        self.inlets.append(inlet.Inlet(self.domain, polygon1, outward_vector1))
+        self.inlets.append(inlet.Inlet(self.domain, polygon1, ep1, outward_vector1))
 
 
@@ -78,10 +81,13 @@
         # Short hands
         w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
-        h = 0.5*self.height*self.culvert_vector    # Vector of length=height in the
+        h = self.apron*self.culvert_vector    # Vector of length=height in the
                              # direction of the culvert
 
+        gap = (1 + self.enquiry_gap)*h
+
         self.inlet_polygons = []
+        self.inlet_equiry_pts = []
 
-        # Build exchange polygon and enquiry points 0 and 1
+        # Build exchange polygon and enquiry point
         for i in [0, 1]:
             i0 = (2*i-1)
@@ -90,12 +96,14 @@
             p2 = p1 + i0*h
             p3 = p0 + i0*h
+            ep = self.end_points[i] + i0*gap
+
             self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
+            self.inlet_equiry_pts.append(ep)
 
         # Check that enquiry points are outside inlet polygons
         for i in [0,1]:
             polygon = self.inlet_polygons[i]
-            # FIXME (SR) Probably should calculate the area of all the triangles
-            # associated with this polygon, as there is likely to be some
-            # inconsistency between triangles and ploygon
+            ep = self.inlet_equiry_pts[i]
+           
             area = polygon_area(polygon)
             
@@ -103,5 +111,7 @@
             msg += ' has area = %f' % area
             assert area > 0.0, msg
-            
+
+            msg = 'Enquiry point falls inside an exchange polygon.'
+            assert not inside_polygon(ep, polygon), msg
     
     def get_inlets(self):
@@ -123,3 +133,8 @@
     
         return self.height
+
+
+    def get_culvert_apron(self):
+
+        return self.apron
         
Index: /trunk/anuga_core/source/anuga/structures/culvert_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/culvert_operator.py	(revision 7983)
+++ /trunk/anuga_core/source/anuga/structures/culvert_operator.py	(revision 7984)
@@ -22,23 +22,38 @@
                  width,
                  height=None,
+                 apron=None,
+                 enquiry_gap=0.2,
                  verbose=False):
         
         self.domain = domain
         self.domain.set_fractional_step_operator(self)
-        end_points = [end_point0, end_point1]
+        self.end_points = [end_point0, end_point1]
         
         if height is None:
             height = width
 
-        self.width = width
+        if apron is None:
+            apron = width
+
+        self.width  = width
         self.height = height
+        self.apron  = apron
+        self.enquiry_gap = enquiry_gap
+        self.verbose = verbose
        
-        self.culvert = Boyd_box_culvert(self.domain, end_points, self.width, self.height)
-        print self.culvert
+        self.culvert = Boyd_box_culvert(self.domain,
+                                        self.end_points,
+                                        self.width,
+                                        self.height,
+                                        self.apron,
+                                        self.enquiry_gap,
+                                        self.verbose)
+        
         self.routine = self.culvert.routine
-        print self.routine
-        self.inlets = self.culvert.get_inlets()
-   
-        self.print_stats()
+
+        self.inlets  = self.culvert.get_inlets()
+
+        if self.verbose:
+            self.print_stats()
 
 
@@ -48,4 +63,5 @@
         
         Q, barrel_speed, outlet_depth = self.routine()
+
 
         inflow  = self.routine.get_inflow()
@@ -54,34 +70,40 @@
 
         old_inflow_height = inflow.get_average_height()
-		old_inflow_xmom = inflow.get_average_xmom()
-		old_inflow_ymom = inflow.get_average_ymom()
-		
-		if old_inflow_height > 0.0 :
-			Qstar = Q/old_inflow_height/inflow.get_area()
-		else:
-			Qstar = 0.0
-
-		factor = 1.0/(1.0 + Qstar*timestep)
-
-		
-		
-		new_inflow_height = old_inflow_height*factor
-		new_inflow_xmom = old_inflow_xmom*factor
-		new_inflow_ymom = old_inflow_ymom*factor
+	old_inflow_xmom = inflow.get_average_xmom()
+	old_inflow_ymom = inflow.get_average_ymom()
+		
+	if old_inflow_height > 0.0 :
+            Qstar = Q/old_inflow_height
+	else:
+            Qstar = 0.0
+
+	factor = 1.0/(1.0 + Qstar*timestep/inflow.get_area())
+
+		
+		
+	new_inflow_height = old_inflow_height*factor
+	new_inflow_xmom = old_inflow_xmom*factor
+	new_inflow_ymom = old_inflow_ymom*factor
 		
 
         inflow.set_heights(new_inflow_height)
+
+        #inflow.set_xmoms(Q/inflow.get_area())
+        #inflow.set_ymoms(0.0)
+
+
         inflow.set_xmoms(new_inflow_xmom)
         inflow.set_ymoms(new_inflow_ymom)
 
-		
-		# set outflow
-		if old_inflow_height > 0.0 :
-			timestep_star = timestep*new_inflow_height/old_inflow_height
-		else:
-			timestep_star = 0.0
-		
-		print Q, barrel_speed, outlet_depth, Qstar, factor, timestep_star
-		
+
+        loss = (old_inflow_height - new_inflow_height)*inflow.get_area()
+
+		
+	# set outflow
+	if old_inflow_height > 0.0 :
+            timestep_star = timestep*new_inflow_height/old_inflow_height
+	else:
+            timestep_star = 0.0
+
 		
         outflow_extra_height = Q*timestep_star/outflow.get_area()
@@ -90,7 +112,23 @@
 		
 
-        outflow.set_heights(outflow.get_average_height() + outflow_extra_height)
-        outflow.set_xmoms(outflow.get_average_xmom() + outflow_extra_momentum[0] )
-        outflow.set_ymoms(outflow.get_average_ymom() + outflow_extra_momentum[1] )
+        gain = outflow_extra_height*outflow.get_area()
+        
+        #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
+        #print '  ', loss, gain
+
+
+        new_outflow_height = outflow.get_average_height() + outflow_extra_height
+        new_outflow_xmom = outflow.get_average_xmom() + outflow_extra_momentum[0]
+        new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1]
+
+        outflow.set_heights(new_outflow_height)
+
+        outflow.set_xmoms(barrel_speed*new_outflow_height*outflow_direction[0])
+        outflow.set_ymoms(barrel_speed*new_outflow_height*outflow_direction[1])
+
+        #outflow.set_xmoms(new_outflow_xmom)
+        #outflow.set_ymoms(new_outflow_ymom)
+        
+        #print '   outflow volume ',outflow.get_total_water_volume()
 
     def print_stats(self):
@@ -99,4 +137,10 @@
         print 'Generic Culvert Operator'
         print '====================================='
+
+        print 'Culvert'
+        print self.culvert
+
+        print 'Culvert Routine'
+        print self.routine
         
         for i, inlet in enumerate(self.inlets):
Index: /trunk/anuga_core/source/anuga/structures/culvert_routine.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/culvert_routine.py	(revision 7983)
+++ /trunk/anuga_core/source/anuga/structures/culvert_routine.py	(revision 7984)
@@ -50,5 +50,5 @@
         # Determine flow direction based on total energy difference
 
-        self.delta_total_energy = self.inlets[0].get_average_total_energy() - self.inlets[1].get_average_total_energy()
+        self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
 
         self.inflow  = self.inlets[0]
Index: /trunk/anuga_core/source/anuga/structures/inlet.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/inlet.py	(revision 7983)
+++ /trunk/anuga_core/source/anuga/structures/inlet.py	(revision 7984)
@@ -9,18 +9,20 @@
     """
 
-    def __init__(self, domain, polygon, outward_culvert_vector=None):
+    def __init__(self, domain, polygon, enquiry_pt,  outward_culvert_vector=None, verbose=False):
 
         self.domain = domain
         self.domain_bounding_polygon = self.domain.get_boundary_polygon()
         self.polygon = polygon
+        self.enquiry_pt = enquiry_pt
         self.outward_culvert_vector = outward_culvert_vector
+        self.verbose = verbose
 
         # FIXME (SR) Using get_triangle_containing_point which needs to be sped up
 
-        self.compute_triangle_indices()
+        self.compute_indices()
         self.compute_area()
 
 
-    def compute_triangle_indices(self):
+    def compute_indices(self):
 
         # Get boundary (in absolute coordinates)
@@ -28,21 +30,20 @@
         domain_centroids = self.domain.get_centroid_coordinates(absolute=True)
 
-        self.inlet_polygon = self.polygon
 
         # Check that polygon lies within the mesh.
-        for point in self.inlet_polygon:
-                msg = 'Point %s in polygon for forcing term' %  str(point)
+        for point in self.polygon + self.enquiry_pt:
+                msg = 'Point %s ' %  str(point)
                 msg += ' did not fall within the domain boundary.'
                 assert is_inside_polygon(point, bounding_polygon), msg
 
 
-        self.triangle_indices = inside_polygon(domain_centroids, self.inlet_polygon, verbose=True)
+        self.triangle_indices = inside_polygon(domain_centroids, self.polygon, verbose=self.verbose)
 
         if len(self.triangle_indices) == 0:
-            region = 'Inlet polygon=%s' % (self.inlet_polygon)
-            msg = 'No triangles have been identified in '
-            msg += 'specified region: %s' % region
+            region = 'Inlet polygon=%s' % (self.polygon)
+            msg = 'No triangles have been identified in region '
             raise Exception, msg
 
+        self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_pt)
 
 
@@ -54,6 +55,5 @@
         if len(self.triangle_indices) == 0:
             region = 'Inlet polygon=%s' % (self.inlet_polygon)
-            msg = 'No triangles have been identified in '
-            msg += 'specified region: %s' % region
+            msg = 'No triangles have been identified in region '
             raise Exception, msg
         
@@ -115,4 +115,5 @@
         return num.sum(self.get_ymoms()*self.get_areas())/self.area
     
+
     def get_heights(self):
     
@@ -175,4 +176,69 @@
 
 
+    def get_enquiry_stage(self):
+
+        return self.domain.quantities['stage'].centroid_values[self.enquiry_index]
+
+
+    def get_enquiry_xmom(self):
+
+        return self.domain.quantities['xmomentum'].centroid_values[self.enquiry_index]
+
+    def get_enquiry_ymom(self):
+
+        return self.domain.quantities['ymomentum'].centroid_values[self.enquiry_index]
+
+
+    def get_enquiry_elevation(self):
+
+        return self.domain.quantities['elevation'].centroid_values[self.enquiry_index]
+
+    def get_enquiry_height(self):
+
+        return self.get_enquiry_stage() - self.get_enquiry_elevation()
+
+
+    def get_enquiry_velocity(self):
+
+            height = self.get_enquiry_height()
+            u = self.get_enquiry_xmom()/(height + velocity_protection/height)
+            v = self.get_enquiry_ymom()/(height + velocity_protection/height)
+
+            return u, v
+
+
+    def get_enquiry_xvelocity(self):
+
+            height = self.get_enquiry_height()
+            return self.get_enquiry_xmom()/(height + velocity_protection/height)
+
+    def get_enquiry_yvelocity(self):
+
+            height = self.get_enquiry_height()
+            return self.get_enquiry_ymom()/(height + velocity_protection/height)
+
+
+    def get_enquiry_speed(self):
+
+            u, v = self.get_enquiry_velocity()
+
+            return math.sqrt(u**2 + v**2)
+
+
+    def get_enquiry_velocity_head(self):
+
+        return 0.5*self.get_enquiry_speed()**2/g
+
+
+    def get_enquiry_total_energy(self):
+
+        return self.get_enquiry_velocity_head() + self.get_enquiry_stage()
+
+
+    def get_enquiry_specific_energy(self):
+
+        return self.get_enquiry_velocity_head() + self.get_enquiry_height()
+
+
     def set_heights(self,height):
 
@@ -187,10 +253,10 @@
     def set_xmoms(self,xmom):
 
-        self.xmoms=self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom)
+        self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom)
 
 
     def set_ymoms(self,ymom):
 
-        self.xmoms=self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom)
+        self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom)
 
 
Index: /trunk/anuga_core/source/anuga/structures/testing_wide_bridge.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/testing_wide_bridge.py	(revision 7983)
+++ /trunk/anuga_core/source/anuga/structures/testing_wide_bridge.py	(revision 7984)
@@ -46,5 +46,5 @@
 
 points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
-                                               len1=length, len2=width)
+                                                    len1=length, len2=width)
 domain = Domain(points, vertices, boundary)   
 domain.set_name('Test_WIDE_BRIDGE')                 # Output name
@@ -168,9 +168,28 @@
 """
 from anuga.structures.culvert_operator import Culvert_operator
+#culvert0 = Culvert_operator(domain,
+#                            end_point0=[40.0, 75.0],
+#                            end_point1=[50.0, 75.0],
+#                            width=50.0,
+#                            height=10.0,
+#                            apron=5.0,
+#                            verbose=False)
+
+
 culvert1 = Culvert_operator(domain,
-                            end_point0=[40.0, 75.0],
-                            end_point1=[50.0, 75.0],
-                            width=50.0,
+                            end_point0=[40.0, 87.5],
+                            end_point1=[50.0, 87.5],
+                            width=25.0,
                             height=10.0,
+                            apron=5.0,
+                            verbose=False)
+
+
+culvert2 = Culvert_operator(domain,
+                            end_point0=[40.0, 62.5],
+                            end_point1=[50.0, 62.5],
+                            width=25.0,
+                            height=10.0,
+                            apron=5.0,
                             verbose=False)
 
@@ -220,4 +239,5 @@
 for t in domain.evolve(yieldstep = 1, finaltime = 100):
     print domain.timestepping_statistics()
+    print domain.volumetric_balance_statistics()
     
 
