Changeset 5723

Timestamp:

Sep 3, 2008, 2:33:25 PM (16 years ago)

Author:

ole

Message:

Refactored get_flow_through_cross_section in preparation for
get_energy_through_cross_section. These function work only on sww files.
Runtime versions still need to be done.

Location:

anuga_core/source/anuga

Files:

• fit_interpolate/interpolate.py (modified) (3 diffs)
• shallow_water/data_manager.py (modified) (8 diffs)

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -721 +721 @@
                     if verbose is True:
                         import sys
-                        if sys.platform == 'win32':
+                        if sys.platform == 'win32': # FIXME (Ole): Why only Windoze?
                             from anuga.utilities.polygon import plot_polygons
                             #out_interp_pts = take(interpolation_points,[indices])
                             title = 'Interpolation points fall outside specified mesh'
-                            plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts],
-                                          ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose)
+                            plot_polygons([mesh_boundary_polygon,
+                                           interpolation_points,
+                                           out_interp_pts],
+                                          ['line','point','outside'],
+                                          figname='points_boundary_out',
+                                          label=title,
+                                          verbose=verbose)
 
                     # Joaquim Luis suggested this as an Exception, so
@@ -735 +740 @@
                         print msg
                     #raise Exception(msg)
+                    
             elif triangles is None and vertex_coordinates is not None:#jj
                 #Dealing with sts file
@@ -748 +754 @@
                     from anuga.utilities.polygon import plot_polygons
                     title = 'Interpolation function: Polygon and interpolation points'
-                    plot_polygons([mesh_boundary_polygon,interpolation_points],
-                                         ['line','point'],figname='points_boundary',label=title,verbose=verbose)
+                    plot_polygons([mesh_boundary_polygon,
+                                   interpolation_points],
+                                  ['line','point'],
+                                  figname='points_boundary',
+                                  label=title,
+                                  verbose=verbose)
 
             m = len(self.interpolation_points)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -6252 +6252 @@
                                       quantities=None,
                                       verbose=False):
-    """Get and rebuild mesh structure and the associated quantities from sww file
-
+    """Get and rebuild mesh structure and associated quantities from sww file
+    
     Input:
         filename - Name os sww file
@@ -6263 +6263 @@
         quantities - arrays with quantity values at each mesh node
         time - vector of stored timesteps
+        
+        
+    This function is used by e.g. get_interpolated_quantities_at_polyline_midpoints
     """
  
@@ -6327 +6330 @@
 
 
-
-def get_flow_through_cross_section(filename,
-                                   polyline,
-                                   verbose=False):
-    """Obtain flow (m^3/s) perpendicular to specified cross section.
-
-    Inputs:
-        filename: Name of sww file
-        polyline: Representation of desired cross section - it may contain multiple
-                  sections allowing for complex shapes. Assume absolute UTM coordinates.
-                  Format [[x0, y0], [x1, y1], ...]
+    
+
+def get_interpolated_quantities_at_polyline_midpoints(filename,
+                                                      quantity_names=None,
+                                                      polyline=None, 
+                                                      verbose=False):
+    """Get values for quantities interpolated to polyline midpoints from sww file
+    
+    Input:
+        filename - Name of sww file
+        quantity_names - Names of quantities to load
+        polyline: Representation of desired cross section - it may contain 
+                  multiple sections allowing for complex shapes. Assume 
+                  absolute UTM coordinates.
+                  Format [[x0, y0], [x1, y1], ...]              
 
     Output:
-        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
-    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.
-
-    """
-
+        segments: list of instances of class Triangle_intersection 
+        interpolation_function: Instance of class Interpolation_function
+          
+          
+      This function is used by get_flow_through_cross_section and 
+      get_energy_through_cross_section          
+    """
+    
     from anuga.fit_interpolate.interpolate import Interpolation_function
-
-    quantity_names =['elevation',
-                     'stage',
-                     'xmomentum',
-                     'ymomentum']
 
     # Get mesh and quantities from sww file
@@ -6371 +6369 @@
     # Find all intersections and associated triangles.
     segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
-    #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 as interpolation_points.
@@ -6388 +6383 @@
         print 'total number of interpolation points = %d'\
               %len(interpolation_points)
-
+              
+              
+    
     I = Interpolation_function(time,
                                quantities,
@@ -6396 +6393 @@
                                interpolation_points=interpolation_points,
                                verbose=verbose)
+                               
+    return segments, I
+    
+
+def get_flow_through_cross_section(filename,
+                                   polyline,
+                                   verbose=False):
+    """Obtain flow (m^3/s) perpendicular to specified cross section.
+
+    Inputs:
+        filename: Name of sww file
+        polyline: Representation of desired cross section - it may contain 
+                  multiple sections allowing for complex shapes. Assume 
+                  absolute UTM coordinates.
+                  Format [[x0, y0], [x1, y1], ...]
+
+    Output:
+        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 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.
+
+    """
+
+    quantity_names =['elevation',
+                     'stage',
+                     'xmomentum',
+                     'ymomentum']
+
+
+    # Get values for quantities at each midpoint of poly line from sww file
+    X = get_interpolated_quantities_at_polyline_midpoints(filename,
+                                                          quantity_names=quantity_names,
+                                                          polyline=polyline,
+                                                          verbose=verbose)    
+    segments, interpolation_function = X
+
+    
+    # Get vectors for time and interpolation_points
+    time = interpolation_function.time
+    interpolation_points = interpolation_function.interpolation_points    
 
     if verbose: print 'Computing hydrograph'
@@ -6403 +6446 @@
         total_flow=0
         for i, p in enumerate(interpolation_points):
-            elevation, stage, uh, vh = I(t, point_id=i)
+            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
             normal = segments[i].normal
 
@@ -6418 +6461 @@
         # Store flow at this timestep    
         Q.append(total_flow)
-