Changeset 5288 for anuga_core/source/anuga/shallow_water/data_manager.py

Timestamp:

May 7, 2008, 4:01:52 PM (16 years ago)

Author:

ole

Message:

Implemeted and tested ticket:175 - flow_through_cross_section.
Also updated manual

File:

• anuga_core/source/anuga/shallow_water/data_manager.py (modified) (1 diff)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -5621 +5621 @@
         filename: Name of sww file
         polyline: Representation of desired cross section - it may contain multiple
-                  sections allowing for complex shapes.
+                  sections allowing for complex shapes. Assume absolute UTM coordinates.
+                  Format [[x0, y0], [x1, y1], ...]
 
     Output:
-        time: All stored times
+        time: All stored times in sww file
         Q: Hydrograph of total flow across given segments for all stored times.
 
     The normal flow is computed for each triangle intersected by the polyline and
-    added up.  multiple segments at different angles are specified the normal flows
+    added up.  Multiple segments at different angles are specified the normal flows
     may partially cancel each other.
 
-    """
+    The typical usage of this function would be to get flow through a channel,
+    and the polyline would then be a cross section perpendicular to the flow.
+
+    """
+
+    from anuga.fit_interpolate.interpolate import Interpolation_function
+
+    quantity_names =['elevation',
+                     'stage',
+                     'xmomentum',
+                     'ymomentum']
 
     # Get mesh and quantities from sww file
     X = get_mesh_and_quantities_from_file(filename,
-                                          quantities=['elevation',
-                                                      'stage',
-                                                      'xmomentum',
-                                                      'ymomentum'], 
+                                          quantities=quantity_names,
                                           verbose=verbose)
     mesh, quantities, time = X
 
 
+    # Adjust polyline to mesh spatial origin
+    polyline = mesh.geo_reference.get_relative(polyline)
     
     # Find all intersections and associated triangles.
-    
-    mesh.get_intersecting_segments(polyline)
+    segments = mesh.get_intersecting_segments(polyline)
+    #print
+    #for x in segments:
+    #    print x
+
     
     # Then store for each triangle the length of the intersecting segment(s),
-    # right hand normal(s) and midpoints.
+    # right hand normal(s) and midpoints as interpolation_points.
+    interpolation_points = []
+    for segment in segments:
+        midpoint = sum(array(segment.segment))/2
+        interpolation_points.append(midpoint)
+
+
+    I = Interpolation_function(time,
+                               quantities,
+                               quantity_names=quantity_names,
+                               vertex_coordinates=mesh.nodes,
+                               triangles=mesh.triangles,
+                               interpolation_points=interpolation_points,
+                               verbose=verbose)
+
+    # Compute hydrograph
+    Q = []
+    for t in time:
+        total_flow=0
+        for i, p in enumerate(interpolation_points):
+            elevation, stage, uh, vh = I(t, point_id=i)
+            normal = segments[i].normal
+
+            # Inner product of momentum vector with segment normal [m^2/s] 
+            normal_momentum = uh*normal[0] + vh*normal[1] 
+
+            # Flow across this segment [m^3/s]
+            segment_flow = normal_momentum*segments[i].length
+
+            # Accumulate
+            total_flow += segment_flow
+             
+
+        # Store flow at this timestep    
+        Q.append(total_flow)
+    
+
 
     return time, Q