Changeset 7770

Timestamp:

Jun 2, 2010, 5:17:47 PM (15 years ago)

Author:

James Hudson

Message:

Broke up data_manager into more modules - some failing tests.

Location:

trunk/anuga_core

Files:

• documentation/user_manual/demos/channel1.py (modified) (3 diffs)
• source/anuga/file/sww.py (modified) (3 diffs)
• source/anuga/file/test_sww.py (modified) (2 diffs)
• source/anuga/file_conversion/file_conversion.py (modified) (1 diff)
• source/anuga/file_conversion/test_file_conversion.py (modified) (2 diffs)
• source/anuga/file_conversion/test_urs2sts.py (added)
• source/anuga/file_conversion/test_urs2sww.py (added)
• source/anuga/file_conversion/urs2sww.py (modified) (4 diffs)
• source/anuga/file_conversion/urs2txt.py (modified) (2 diffs)
• source/anuga/pmesh/mesh_interface.py (modified) (5 diffs)
• source/anuga/shallow_water/data_manager.py (modified) (2 diffs)
• source/anuga/shallow_water/sww_interrogate.py (added)
• source/anuga/shallow_water/test_data_manager.py (modified) (2 diffs)
• source/anuga/shallow_water/test_shallow_water_domain.py (modified) (3 diffs)

trunk/anuga_core/documentation/user_manual/demos/channel1.py

@@ -7 +7 @@
 # Import necessary modules
 #------------------------------------------------------------------------------
+
+# Import standard shallow water domain and standard boundaries.
+import anuga
+
+# Import rectangular cross mesh generator
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-from anuga.shallow_water import Domain
-from anuga.shallow_water import Reflective_boundary
-from anuga.shallow_water import Dirichlet_boundary
+
 
 #------------------------------------------------------------------------------
@@ -18 +21 @@
                                                len1=10.0, len2=5.0) # Mesh
 
-domain = Domain(points, vertices, boundary)  # Create domain
+domain = anuga.Domain(points, vertices, boundary)  # Create domain
 domain.set_name('channel1')                  # Output name
 
@@ -35 +38 @@
 # Setup boundary conditions
 #------------------------------------------------------------------------------
-Bi = Dirichlet_boundary([0.4, 0, 0])         # Inflow
-Br = Reflective_boundary(domain)             # Solid reflective wall
+Bi = anuga.Dirichlet_boundary([0.4, 0, 0])         # Inflow
+Br = anuga.Reflective_boundary(domain)             # Solid reflective wall
 
 domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})

trunk/anuga_core/source/anuga/file/sww.py

@@ -1196 +1196 @@
 
 ##
+# @brief Get mesh and quantity data from an SWW file.
+# @param filename Path to data file to read.
+# @param quantities UNUSED!
+# @param verbose True if this function is to be verbose.
+# @return (mesh, quantities, time) where mesh is the mesh data, quantities is
+#         a dictionary of {name: value}, and time is the time vector.
+# @note Quantities extracted: 'elevation', 'stage', 'xmomentum' and 'ymomentum'
+def get_mesh_and_quantities_from_file(filename,
+                                      quantities=None,
+                                      verbose=False):
+    """Get and rebuild mesh structure and associated quantities from sww file
+
+    Input:
+        filename - Name os sww file
+        quantities - Names of quantities to load
+
+    Output:
+        mesh - instance of class Interpolate
+               (including mesh and interpolation functionality)
+        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
+    """
+
+    # FIXME (Ole): Maybe refactor filefunction using this more fundamental code.
+
+    import types
+    from Scientific.IO.NetCDF import NetCDFFile
+    from shallow_water import Domain
+    from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
+
+    if verbose: log.critical('Reading from %s' % filename)
+
+    fid = NetCDFFile(filename, netcdf_mode_r)    # Open existing file for read
+    time = fid.variables['time'][:]    # Time vector
+    time += fid.starttime[0]
+
+    # Get the variables as numeric arrays
+    x = fid.variables['x'][:]                   # x-coordinates of nodes
+    y = fid.variables['y'][:]                   # y-coordinates of nodes
+    elevation = fid.variables['elevation'][:]   # Elevation
+    stage = fid.variables['stage'][:]           # Water level
+    xmomentum = fid.variables['xmomentum'][:]   # Momentum in the x-direction
+    ymomentum = fid.variables['ymomentum'][:]   # Momentum in the y-direction
+
+    # Mesh (nodes (Mx2), triangles (Nx3))
+    nodes = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]), axis=1)
+    triangles = fid.variables['volumes'][:]
+
+    # Get geo_reference
+    try:
+        geo_reference = Geo_reference(NetCDFObject=fid)
+    except: #AttributeError, e:
+        # Sww files don't have to have a geo_ref
+        geo_reference = None
+
+    if verbose: log.critical('    building mesh from sww file %s' % filename)
+
+    boundary = None
+
+    #FIXME (Peter Row): Should this be in mesh?
+    if fid.smoothing != 'Yes':
+        nodes = nodes.tolist()
+        triangles = triangles.tolist()
+        nodes, triangles, boundary = weed(nodes, triangles, boundary)
+
+    try:
+        mesh = Mesh(nodes, triangles, boundary, geo_reference=geo_reference)
+    except AssertionError, e:
+        fid.close()
+        msg = 'Domain could not be created: %s. "' % e
+        raise DataDomainError, msg
+
+    quantities = {}
+    quantities['elevation'] = elevation
+    quantities['stage'] = stage
+    quantities['xmomentum'] = xmomentum
+    quantities['ymomentum'] = ymomentum
+
+    fid.close()
+
+    return mesh, quantities, time
+
+
+
+##
 # @brief 
 # @parm time 
@@ -1251 +1339 @@
 
 
-
-
 ##
 # @brief 
@@ -1310 +1396 @@
 
     return coordinates, volumes, boundary
+
+
+

trunk/anuga_core/source/anuga/file/test_sww.py

@@ -8 +8 @@
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
 from anuga.shallow_water.shallow_water_domain import Domain
-from sww import load_sww_as_domain
+from sww import load_sww_as_domain, weed, get_mesh_and_quantities_from_file
 
 # boundary functions
@@ -158 +158 @@
 
 
+    def test_get_mesh_and_quantities_from_sww_file(self):
+        """test_get_mesh_and_quantities_from_sww_file(self):
+        """     
+        
+        # Generate a test sww file with non trivial georeference
+        
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        # Setup
+        from mesh_factory import rectangular
+
+        # Create basic mesh (100m x 5m)
+        width = 5
+        length = 50
+        t_end = 10
+        points, vertices, boundary = rectangular(length, width, 50, 5)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary,
+                        geo_reference = Geo_reference(56,308500,6189000))
+
+        domain.set_name('flowtest')
+        swwfile = domain.get_name() + '.sww'
+        domain.set_datadir('.')
+
+        Br = Reflective_boundary(domain)    # Side walls
+        Bd = Dirichlet_boundary([1, 0, 0])  # inflow
+
+        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
+
+        for t in domain.evolve(yieldstep=1, finaltime = t_end):
+            pass
+
+        
+        # Read it
+
+        # Get mesh and quantities from sww file
+        X = get_mesh_and_quantities_from_file(swwfile,
+                                              quantities=['elevation',
+                                                          'stage',
+                                                          'xmomentum',
+                                                          'ymomentum'], 
+                                              verbose=False)
+        mesh, quantities, time = X
+        
+
+        # Check that mesh has been recovered
+        assert num.alltrue(mesh.triangles == domain.get_triangles())
+        assert num.allclose(mesh.nodes, domain.get_nodes())
+
+        # Check that time has been recovered
+        assert num.allclose(time, range(t_end+1))
+
+        # Check that quantities have been recovered
+        # (sww files use single precision)
+        z=domain.get_quantity('elevation').get_values(location='unique vertices')
+        assert num.allclose(quantities['elevation'], z)
+
+        for q in ['stage', 'xmomentum', 'ymomentum']:
+            # Get quantity at last timestep
+            q_ref=domain.get_quantity(q).get_values(location='unique vertices')
+
+            #print q,quantities[q]
+            q_sww=quantities[q][-1,:]
+
+            msg = 'Quantity %s failed to be recovered' %q
+            assert num.allclose(q_ref, q_sww, atol=1.0e-6), msg
+            
+        # Cleanup
+        os.remove(swwfile)
+        
+        
+
+    def test_weed(self):
+        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
+        volumes1 = [[0,1,2],[3,4,5]]
+        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
+        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
+
+        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
+
+        assert len(points2)==len(coordinates2)
+        for i in range(len(coordinates2)):
+            coordinate = tuple(coordinates2[i])
+            assert points2.has_key(coordinate)
+            points2[coordinate]=i
+
+        for triangle in volumes1:
+            for coordinate in triangle:
+                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
+                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
+
 #################################################################################
 

trunk/anuga_core/source/anuga/file_conversion/file_conversion.py

@@ -488 +488 @@
 
 
+##
+# @brief 
+# @param filename 
+# @param x 
+# @param y 
+# @param z 
+def write_obj(filename, x, y, z):
+    """Store x,y,z vectors into filename (obj format).
+
+       Vectors are assumed to have dimension (M,3) where
+       M corresponds to the number elements.
+       triangles are assumed to be disconnected
+
+       The three numbers in each vector correspond to three vertices,
+
+       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
+    """
+
+    import os.path
+
+    root, ext = os.path.splitext(filename)
+    if ext == '.obj':
+        FN = filename
+    else:
+        FN = filename + '.obj'
+
+    outfile = open(FN, 'wb')
+    outfile.write("# Triangulation as an obj file\n")
+
+    M, N = x.shape
+    assert N == 3  #Assuming three vertices per element
+
+    for i in range(M):
+        for j in range(N):
+            outfile.write("v %f %f %f\n" % (x[i,j], y[i,j], z[i,j]))
+
+    for i in range(M):
+        base = i * N
+        outfile.write("f %d %d %d\n" % (base+1, base+2, base+3))
+
+    outfile.close()
+

trunk/anuga_core/source/anuga/file_conversion/test_file_conversion.py

@@ -8 +8 @@
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
 from anuga.shallow_water.sww_file import SWW_file
-from data_manager import extent_sww
+from anuga.file.sww import extent_sww
 from anuga.config import netcdf_float, epsilon, g
 from Scientific.IO.NetCDF import NetCDFFile
 from anuga.file_conversion.file_conversion import tsh2sww, \
                         pmesh_to_domain_instance
+
+from anuga.file.mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
+                            NORTH_VELOCITY_LABEL
 
 import sys
@@ -1570 +1573 @@
         #print "sww_file",sww_file
         #dem_file = tempfile.mktemp(".dem")
-        domain = sww2domain(sww_file) ###, dem_file)
+        domain = load_sww_as_domain(sww_file) ###, dem_file)
         domain.check_integrity()
 

trunk/anuga_core/source/anuga/file_conversion/urs2sww.py

@@ +1 @@
+import numpy as num
+from Scientific.IO.NetCDF import NetCDFFile
+
+from anuga.file.urs import Read_urs
+
+from anuga.coordinate_transforms.redfearn import redfearn, \
+     convert_from_latlon_to_utm
+
+from anuga.geospatial_data.geospatial_data import ensure_absolute, \
+                                                    Geospatial_data
+
+from mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
+                            NORTH_VELOCITY_LABEL
+                            
+from anuga.utilities.numerical_tools import ensure_numeric                            
+
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
+                            netcdf_float
+
+from sww_file import Read_sww, Write_sww
+
+
 ################################################################################
 # CONVERTING UNGRIDDED URS DATA TO AN SWW FILE
 ################################################################################
-
-WAVEHEIGHT_MUX_LABEL = '-z-mux'
-EAST_VELOCITY_LABEL =  '-e-mux'
-NORTH_VELOCITY_LABEL =  '-n-mux'
 
 ##
@@ -95 +113 @@
     mux = {}
     for quantity, file in map(None, quantities, files_in):
-        mux[quantity] = Urs_points(file)
+        mux[quantity] = Read_urs(file)
 
     # Could check that the depth is the same. (hashing)
@@ -143 +161 @@
     for i in range(a_mux.time_step_count):
         mux_times.append(a_mux.time_step * i)
-    (mux_times_start_i, mux_times_fin_i) = mux2sww_time(mux_times, mint, maxt)
+    (mux_times_start_i, mux_times_fin_i) = read_time_from_mux(mux_times, mint, maxt)
     times = mux_times[mux_times_start_i:mux_times_fin_i]
 
@@ -209 +227 @@
 
     # Do some conversions while writing the sww file
+
+
+def read_time_from_mux(mux_times, mint, maxt):
+    """
+        Read a list of mux times.
+        Return start and finish times which lie within the passed time period.
+    """
+
+    if mint == None:
+        mux_times_start_i = 0
+    else:
+        mux_times_start_i = num.searchsorted(mux_times, mint)
+
+    if maxt == None:
+        mux_times_fin_i = len(mux_times)
+    else:
+        maxt += 0.5 # so if you specify a time where there is
+                    # data that time will be included
+        mux_times_fin_i = num.searchsorted(mux_times, maxt)
+
+    return mux_times_start_i, mux_times_fin_i
+

trunk/anuga_core/source/anuga/file_conversion/urs2txt.py

@@ -1 @@
-
+from anuga.file.urs import Read_urs
 
 ##
@@ -20 +20 @@
     mux = {}
     for quantity, file in map(None, quantities, files_in):
-        mux[quantity] = Urs_points(file)
+        mux[quantity] = Read_urs(file)
 
     # Could check that the depth is the same. (hashing)

trunk/anuga_core/source/anuga/pmesh/mesh_interface.py

@@ -26 +26 @@
                              interior_regions=None,
                              interior_holes=None,
+                             hole_tags=None,
                              poly_geo_reference=None,
                              mesh_geo_reference=None,
@@ -53 +54 @@
     throw an error
     
-    Interior_holes is a list of ploygons for each hole.
+    Interior_holes is a list of polygons for each hole.
+    hole_tags is an optional list of boundary tags for the holes, see
+                boundary_tags parameter.
 
     This function does not allow segments to share points - use underlying
@@ -92 +95 @@
               'interior_regions': interior_regions,
               'interior_holes': interior_holes,
+              'hole_tags': hole_tags,
               'poly_geo_reference': poly_geo_reference,
               'mesh_geo_reference': mesh_geo_reference,
@@ -127 +131 @@
                               interior_regions=None,
                               interior_holes=None,
+                              hole_tags=None,
                               poly_geo_reference=None,
                               mesh_geo_reference=None,
@@ -311 +316 @@
         for polygon in interior_holes:
             m.add_hole_from_polygon(polygon,
+                                    segment_tags=hole_tags,
                                     geo_reference=poly_geo_reference)
        

trunk/anuga_core/source/anuga/shallow_water/data_manager.py

@@ -268 +268 @@
 
 
-##
-# @brief 
-# @param filename 
-# @param x 
-# @param y 
-# @param z 
-def write_obj(filename, x, y, z):
-    """Store x,y,z vectors into filename (obj format).
-
-       Vectors are assumed to have dimension (M,3) where
-       M corresponds to the number elements.
-       triangles are assumed to be disconnected
-
-       The three numbers in each vector correspond to three vertices,
-
-       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
-    """
-
-    import os.path
-
-    root, ext = os.path.splitext(filename)
-    if ext == '.obj':
-        FN = filename
-    else:
-        FN = filename + '.obj'
-
-    outfile = open(FN, 'wb')
-    outfile.write("# Triangulation as an obj file\n")
-
-    M, N = x.shape
-    assert N == 3  #Assuming three vertices per element
-
-    for i in range(M):
-        for j in range(N):
-            outfile.write("v %f %f %f\n" % (x[i,j], y[i,j], z[i,j]))
-
-    for i in range(M):
-        base = i * N
-        outfile.write("f %d %d %d\n" % (base+1, base+2, base+3))
-
-    outfile.close()
 
 ##
@@ -949 +908 @@
         log.critical(msg)
 
-################################################################################
-# Functions to obtain diagnostics from sww files
-################################################################################
-
-
-
-##
-# @brief Get values for quantities interpolated to polyline midpoints from SWW.
-# @param filename Path to file to read.
-# @param quantity_names Quantity names to get.
-# @param polyline Representation of desired cross-section.
-# @param verbose True if this function is to be verbose.
-# @return (segments, i_func) where segments is a list of Triangle_intersection
-#         instances and i_func is an instance of Interpolation_function.
-# @note For 'polyline' assume absolute UTM coordinates.
-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
-
-    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:
-        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
-
-    # Get mesh and quantities from sww file
-    X = get_mesh_and_quantities_from_file(filename,
-                                          quantities=quantity_names,
-                                          verbose=verbose)
-    mesh, quantities, time = X
-
-    # Find all intersections and associated triangles.
-    segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
-
-    # Get midpoints
-    interpolation_points = segment_midpoints(segments)
-
-    # Interpolate
-    if verbose:
-        log.critical('Interpolating - total number of interpolation points = %d'
-                     % len(interpolation_points))
-
-    I = Interpolation_function(time,
-                               quantities,
-                               quantity_names=quantity_names,
-                               vertex_coordinates=mesh.nodes,
-                               triangles=mesh.triangles,
-                               interpolation_points=interpolation_points,
-                               verbose=verbose)
-
-    return segments, I
-
-
-##
-# @brief Obtain flow (m^3/s) perpendicular to specified cross section.
-# @param filename Path to file to read.
-# @param polyline Representation of desired cross-section.
-# @param verbose Trie if this function is to be verbose.
-# @return (time, Q) where time and Q are lists of time and flow respectively.
-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: log.critical('Computing hydrograph')
-
-    # Compute hydrograph
-    Q = []
-    for t in time:
-        total_flow = 0
-        for i in range(len(interpolation_points)):
-            elevation, stage, uh, vh = interpolation_function(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
-
-
-##
-# @brief Get average energy across a cross-section.
-# @param filename Path to file of interest.
-# @param polyline Representation of desired cross-section.
-# @param kind Select energy to compute: 'specific' or 'total'.
-# @param verbose True if this function is to be verbose.
-# @return (time, E) where time and E are lists of timestep and energy.
-def get_energy_through_cross_section(filename,
-                                     polyline,
-                                     kind='total',
-                                     verbose=False):
-    """Obtain average energy head [m] across specified cross section.
-
-    Inputs:
-        polyline: Representation of desired cross section - it may contain
-                  multiple sections allowing for complex shapes. Assume
-                  absolute UTM coordinates.
-                  Format [[x0, y0], [x1, y1], ...]
-        kind:     Select which energy to compute.
-                  Options are 'specific' and 'total' (default)
-
-    Output:
-        E: Average energy [m] across given segments for all stored times.
-
-    The average velocity is computed for each triangle intersected by
-    the polyline and averaged weighted by segment lengths.
-
-    The typical usage of this function would be to get average energy of
-    flow in a channel, and the polyline would then be a cross section
-    perpendicular to the flow.
-
-    #FIXME (Ole) - need name for this energy reflecting that its dimension
-    is [m].
-    """
-
-    from anuga.config import g, epsilon, velocity_protection as h0
-
-    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: log.critical('Computing %s energy' % kind)
-
-    # Compute total length of polyline for use with weighted averages
-    total_line_length = 0.0
-    for segment in segments:
-        total_line_length += segment.length
-
-    # Compute energy
-    E = []
-    for t in time:
-        average_energy = 0.0
-        for i, p in enumerate(interpolation_points):
-            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
-
-            # Depth
-            h = depth = stage-elevation
-
-            # Average velocity across this segment
-            if h > epsilon:
-                # Use protection against degenerate velocities
-                u = uh / (h + h0/h)
-                v = vh / (h + h0/h)
-            else:
-                u = v = 0.0
-
-            speed_squared = u*u + v*v
-            kinetic_energy = 0.5 * speed_squared / g
-
-            if kind == 'specific':
-                segment_energy = depth + kinetic_energy
-            elif kind == 'total':
-                segment_energy = stage + kinetic_energy
-            else:
-                msg = 'Energy kind must be either "specific" or "total". '
-                msg += 'I got %s' % kind
-
-            # Add to weighted average
-            weigth = segments[i].length / total_line_length
-            average_energy += segment_energy * weigth
-
-        # Store energy at this timestep
-        E.append(average_energy)
-
-    return time, E
-
-
-##
-# @brief Return highest elevation where depth > 0.
-# @param filename Path to SWW file of interest.
-# @param polygon If specified resrict to points inside this polygon.
-# @param time_interval If specified resrict to within the time specified.
-# @param verbose True if this function is  to be verbose.
-def get_maximum_inundation_elevation(filename,
-                                     polygon=None,
-                                     time_interval=None,
-                                     verbose=False):
-    """Return highest elevation where depth > 0
-
-    Usage:
-    max_runup = get_maximum_inundation_elevation(filename,
-                                                 polygon=None,
-                                                 time_interval=None,
-                                                 verbose=False)
-
-    filename is a NetCDF sww file containing ANUGA model output.
-    Optional arguments polygon and time_interval restricts the maximum
-    runup calculation
-    to a points that lie within the specified polygon and time interval.
-
-    If no inundation is found within polygon and time_interval the return value
-    is None signifying "No Runup" or "Everything is dry".
-
-    See general function get_maximum_inundation_data for details.
-    """
-
-    runup, _ = get_maximum_inundation_data(filename,
-                                           polygon=polygon,
-                                           time_interval=time_interval,
-                                           verbose=verbose)
-    return runup
-
-
-##
-# @brief Return location of highest elevation where h > 0
-# @param filename Path to SWW file to read.
-# @param polygon If specified resrict to points inside this polygon.
-# @param time_interval If specified resrict to within the time specified.
-# @param verbose True if this function is  to be verbose.
-def get_maximum_inundation_location(filename,
-                                    polygon=None,
-                                    time_interval=None,
-                                    verbose=False):
-    """Return location of highest elevation where h > 0
-
-    Usage:
-    max_runup_location = get_maximum_inundation_location(filename,
-                                                         polygon=None,
-                                                         time_interval=None,
-                                                         verbose=False)
-
-    filename is a NetCDF sww file containing ANUGA model output.
-    Optional arguments polygon and time_interval restricts the maximum
-    runup calculation
-    to a points that lie within the specified polygon and time interval.
-
-    If no inundation is found within polygon and time_interval the return value
-    is None signifying "No Runup" or "Everything is dry".
-
-    See general function get_maximum_inundation_data for details.
-    """
-
-    _, max_loc = get_maximum_inundation_data(filename,
-                                             polygon=polygon,
-                                             time_interval=time_interval,
-                                             verbose=verbose)
-    return max_loc
-
-
-##
-# @brief Compute maximum run up height from SWW file.
-# @param filename Path to SWW file to read.
-# @param polygon If specified resrict to points inside this polygon.
-# @param time_interval If specified resrict to within the time specified.
-# @param use_centroid_values 
-# @param verbose True if this function is to be verbose.
-# @return (maximal_runup, maximal_runup_location)
-def get_maximum_inundation_data(filename, polygon=None, time_interval=None,
-                                use_centroid_values=False,
-                                verbose=False):
-    """Compute maximum run up height from sww file.
-
-    Usage:
-    runup, location = get_maximum_inundation_data(filename,
-                                                  polygon=None,
-                                                  time_interval=None,
-                                                  verbose=False)
-
-    Algorithm is as in get_maximum_inundation_elevation from
-    shallow_water_domain except that this function works with the sww file and
-    computes the maximal runup height over multiple timesteps.
-
-    Optional arguments polygon and time_interval restricts the maximum runup
-    calculation to a points that lie within the specified polygon and time
-    interval.
-
-    Polygon is assumed to be in (absolute) UTM coordinates in the same zone
-    as domain.
-
-    If no inundation is found within polygon and time_interval the return value
-    is None signifying "No Runup" or "Everything is dry".
-    """
-
-    # We are using nodal values here as that is what is stored in sww files.
-
-    # Water depth below which it is considered to be 0 in the model
-    # FIXME (Ole): Allow this to be specified as a keyword argument as well
-
-    from anuga.geometry.polygon import inside_polygon
-    from anuga.config import minimum_allowed_height
-    from Scientific.IO.NetCDF import NetCDFFile
-
-    dir, base = os.path.split(filename)
-
-    iterate_over = get_all_swwfiles(dir, base)
-
-    # Read sww file
-    if verbose: log.critical('Reading from %s' % filename)
-    # FIXME: Use general swwstats (when done)
-
-    maximal_runup = None
-    maximal_runup_location = None
-
-    for file, swwfile in enumerate (iterate_over):
-        # Read sww file
-        filename = join(dir, swwfile+'.sww')
-
-        if verbose: log.critical('Reading from %s' % filename)
-        # FIXME: Use general swwstats (when done)
-
-        fid = NetCDFFile(filename)
-
-        # Get geo_reference
-        # sww files don't have to have a geo_ref
-        try:
-            geo_reference = Geo_reference(NetCDFObject=fid)
-        except AttributeError, e:
-            geo_reference = Geo_reference() # Default georef object
-
-        xllcorner = geo_reference.get_xllcorner()
-        yllcorner = geo_reference.get_yllcorner()
-        zone = geo_reference.get_zone()
-
-        # Get extent
-        volumes = fid.variables['volumes'][:]
-        x = fid.variables['x'][:] + xllcorner
-        y = fid.variables['y'][:] + yllcorner
-
-        # Get the relevant quantities (Convert from single precison)
-        elevation = num.array(fid.variables['elevation'][:], num.float)
-        stage = num.array(fid.variables['stage'][:], num.float)
-
-        # Here's where one could convert nodal information to centroid
-        # information but is probably something we need to write in C.
-        # Here's a Python thought which is NOT finished!!!
-        if use_centroid_values is True:
-            x = get_centroid_values(x, volumes)
-            y = get_centroid_values(y, volumes)
-            elevation = get_centroid_values(elevation, volumes)
-
-        # Spatial restriction
-        if polygon is not None:
-            msg = 'polygon must be a sequence of points.'
-            assert len(polygon[0]) == 2, msg
-            # FIXME (Ole): Make a generic polygon input check in polygon.py
-            # and call it here
-            points = num.ascontiguousarray(num.concatenate((x[:,num.newaxis],
-                                                            y[:,num.newaxis]),
-                                                            axis=1))
-            point_indices = inside_polygon(points, polygon)
-
-            # Restrict quantities to polygon
-            elevation = num.take(elevation, point_indices, axis=0)
-            stage = num.take(stage, point_indices, axis=1)
-
-            # Get info for location of maximal runup
-            points_in_polygon = num.take(points, point_indices, axis=0)
-
-            x = points_in_polygon[:,0]
-            y = points_in_polygon[:,1]
-        else:
-            # Take all points
-            point_indices = num.arange(len(x))
-
-        # Temporal restriction
-        time = fid.variables['time'][:]
-        all_timeindices = num.arange(len(time))
-        if time_interval is not None:
-            msg = 'time_interval must be a sequence of length 2.'
-            assert len(time_interval) == 2, msg
-            msg = 'time_interval %s must not be decreasing.' % time_interval
-            assert time_interval[1] >= time_interval[0], msg
-            msg = 'Specified time interval [%.8f:%.8f] ' % tuple(time_interval)
-            msg += 'must does not match model time interval: [%.8f, %.8f]\n' \
-                   % (time[0], time[-1])
-            if time_interval[1] < time[0]: raise ValueError(msg)
-            if time_interval[0] > time[-1]: raise ValueError(msg)
-
-            # Take time indices corresponding to interval (& is bitwise AND)
-            timesteps = num.compress((time_interval[0] <= time) \
-                                     & (time <= time_interval[1]),
-                                     all_timeindices)
-
-            msg = 'time_interval %s did not include any model timesteps.' \
-                  % time_interval
-            assert not num.alltrue(timesteps == 0), msg
-        else:
-            # Take them all
-            timesteps = all_timeindices
-
-        fid.close()
-
-        # Compute maximal runup for each timestep
-        #maximal_runup = None
-        #maximal_runup_location = None
-        #maximal_runups = [None]
-        #maximal_runup_locations = [None]
-
-        for i in timesteps:
-            if use_centroid_values is True:
-                stage_i = get_centroid_values(stage[i,:], volumes)
-            else:
-                stage_i = stage[i,:]
-
-            depth = stage_i - elevation
-
-            # Get wet nodes i.e. nodes with depth>0 within given region
-            # and timesteps
-            wet_nodes = num.compress(depth > minimum_allowed_height,
-                                     num.arange(len(depth)))
-
-            if num.alltrue(wet_nodes == 0):
-                runup = None
-            else:
-                # Find maximum elevation among wet nodes
-                wet_elevation = num.take(elevation, wet_nodes, axis=0)
-                runup_index = num.argmax(wet_elevation)
-                runup = max(wet_elevation)
-                assert wet_elevation[runup_index] == runup       # Must be True
-
-            if runup > maximal_runup:
-                maximal_runup = runup      # works even if maximal_runup is None
-
-                # Record location
-                wet_x = num.take(x, wet_nodes, axis=0)
-                wet_y = num.take(y, wet_nodes, axis=0)
-                maximal_runup_location = [wet_x[runup_index],wet_y[runup_index]]
-
-    return maximal_runup, maximal_runup_location
 
 

trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -1351 +1351 @@
         assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
 
-
-
-    def test_weed(self):
-        from data_manager import weed
-
-        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
-        volumes1 = [[0,1,2],[3,4,5]]
-        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
-        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
-
-        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
-
-        assert len(points2)==len(coordinates2)
-        for i in range(len(coordinates2)):
-            coordinate = tuple(coordinates2[i])
-            assert points2.has_key(coordinate)
-            points2[coordinate]=i
-
-        for triangle in volumes1:
-            for coordinate in triangle:
-                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
-                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
 
 
@@ -2488 +2466 @@
         
 
-    def test_get_maximum_inundation(self):
-        """Test that sww information can be converted correctly to maximum
-        runup elevation and location (without and with georeferencing)
-
-        This test creates a slope and a runup which is maximal (~11m) at around 10s
-        and levels out to the boundary condition (1m) at about 30s.
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Setup
-
-        from mesh_factory import rectangular
-
-        # Create basic mesh (100m x 100m)
-        points, vertices, boundary = rectangular(20, 5, 100, 50)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.default_order = 2
-        domain.set_minimum_storable_height(0.01)
-
-        domain.set_name('runuptest')
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-
-        # FIXME (Ole): Backwards compatibility
-        # Look at sww file and see what happens when
-        # domain.tight_slope_limiters = 1
-        domain.tight_slope_limiters = 0
-        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)        
-        
-        Br = Reflective_boundary(domain)
-        Bd = Dirichlet_boundary([1.0,0,0])
-
-
-        #---------- First run without geo referencing
-        
-        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
-        domain.set_quantity('stage', -6)
-        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = 50):
-            pass
-
-
-        # Check maximal runup
-        runup = get_maximum_inundation_elevation(swwfile)
-        location = get_maximum_inundation_location(swwfile)
-        #print 'Runup, location', runup, location
-        assert num.allclose(runup, 11) or num.allclose(runup, 12) # old limiters
-        assert num.allclose(location[0], 15) or num.allclose(location[0], 10)
-
-        # Check final runup
-        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
-        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
-        # print 'Runup, location:',runup, location        
-        assert num.allclose(runup, 1)
-        assert num.allclose(location[0], 65)
-
-        # Check runup restricted to a polygon
-        p = [[50,1], [99,1], [99,49], [50,49]]
-        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
-        location = get_maximum_inundation_location(swwfile, polygon=p)
-        #print runup, location        
-        assert num.allclose(runup, 4)
-        assert num.allclose(location[0], 50)                
-
-        # Check that mimimum_storable_height works
-        fid = NetCDFFile(swwfile, netcdf_mode_r) # Open existing file
-        
-        stage = fid.variables['stage'][:]
-        z = fid.variables['elevation'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]        
-
-        
-        
-        for i in range(stage.shape[0]):
-            h = stage[i]-z # depth vector at time step i
-            
-            # Check every node location
-            for j in range(stage.shape[1]):
-                # Depth being either exactly zero implies
-                # momentum being zero.
-                # Or else depth must be greater than or equal to
-                # the minimal storable height
-                if h[j] == 0.0:
-                    assert xmomentum[i,j] == 0.0
-                    assert ymomentum[i,j] == 0.0                
-                else:
-                    assert h[j] >= domain.minimum_storable_height
-        
-        fid.close()
-
-        # Cleanup
-        os.remove(swwfile)
-        
-
-
-        #------------- Now the same with georeferencing
-
-        domain.time=0.0
-        E = 308500
-        N = 6189000
-        #E = N = 0
-        domain.geo_reference = Geo_reference(56, E, N)
-
-        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
-        domain.set_quantity('stage', -6)
-        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = 50):
-            pass
-
-        # Check maximal runup
-        runup = get_maximum_inundation_elevation(swwfile)
-        location = get_maximum_inundation_location(swwfile)
-        assert num.allclose(runup, 11) or num.allclose(runup, 12) # old limiters
-        assert num.allclose(location[0], 15+E) or num.allclose(location[0], 10+E)
-
-        # Check final runup
-        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
-        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
-        assert num.allclose(runup, 1)
-        assert num.allclose(location[0], 65+E)
-
-        # Check runup restricted to a polygon
-        p = num.array([[50,1], [99,1], [99,49], [50,49]], num.int) + num.array([E, N], num.int)      #array default#
-
-        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
-        location = get_maximum_inundation_location(swwfile, polygon=p)
-        assert num.allclose(runup, 4)
-        assert num.allclose(location[0], 50+E)                
-
-
-        # Cleanup
-        os.remove(swwfile)
-
-
-    def test_get_mesh_and_quantities_from_sww_file(self):
-        """test_get_mesh_and_quantities_from_sww_file(self):
-        """     
-        
-        # Generate a test sww file with non trivial georeference
-        
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        # Setup
-        from mesh_factory import rectangular
-
-        # Create basic mesh (100m x 5m)
-        width = 5
-        length = 50
-        t_end = 10
-        points, vertices, boundary = rectangular(length, width, 50, 5)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary,
-                        geo_reference = Geo_reference(56,308500,6189000))
-
-        domain.set_name('flowtest')
-        swwfile = domain.get_name() + '.sww'
-        domain.set_datadir('.')
-
-        Br = Reflective_boundary(domain)    # Side walls
-        Bd = Dirichlet_boundary([1, 0, 0])  # inflow
-
-        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = t_end):
-            pass
-
-        
-        # Read it
-
-        # Get mesh and quantities from sww file
-        X = get_mesh_and_quantities_from_file(swwfile,
-                                              quantities=['elevation',
-                                                          'stage',
-                                                          'xmomentum',
-                                                          'ymomentum'], 
-                                              verbose=False)
-        mesh, quantities, time = X
-        
-
-        # Check that mesh has been recovered
-        assert num.alltrue(mesh.triangles == domain.get_triangles())
-        assert num.allclose(mesh.nodes, domain.get_nodes())
-
-        # Check that time has been recovered
-        assert num.allclose(time, range(t_end+1))
-
-        # Check that quantities have been recovered
-        # (sww files use single precision)
-        z=domain.get_quantity('elevation').get_values(location='unique vertices')
-        assert num.allclose(quantities['elevation'], z)
-
-        for q in ['stage', 'xmomentum', 'ymomentum']:
-            # Get quantity at last timestep
-            q_ref=domain.get_quantity(q).get_values(location='unique vertices')
-
-            #print q,quantities[q]
-            q_sww=quantities[q][-1,:]
-
-            msg = 'Quantity %s failed to be recovered' %q
-            assert num.allclose(q_ref, q_sww, atol=1.0e-6), msg
-            
-        # Cleanup
-        os.remove(swwfile)
-        
-
-    def test_get_flow_through_cross_section(self):
-        """test_get_flow_through_cross_section(self):
-
-        Test that the total flow through a cross section can be
-        correctly obtained from an sww file.
-        
-        This test creates a flat bed with a known flow through it and tests
-        that the function correctly returns the expected flow.
-
-        The specifics are
-        u = 2 m/s
-        h = 1 m
-        w = 3 m (width of channel)
-
-        q = u*h*w = 6 m^3/s
-
-        #---------- First run without geo referencing        
-        
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        # Setup
-        from mesh_factory import rectangular
-
-        # Create basic mesh (20m x 3m)
-        width = 3
-        length = 20
-        t_end = 3
-        points, vertices, boundary = rectangular(length, width,
-                                                 length, width)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.default_order = 2
-        domain.set_minimum_storable_height(0.01)
-
-        domain.set_name('flowtest')
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-
-        h = 1.0
-        u = 2.0
-        uh = u*h
-
-        Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 3 m inlet: 
-
-
-        
-        domain.set_quantity('elevation', 0.0)
-        domain.set_quantity('stage', h)
-        domain.set_quantity('xmomentum', uh)
-        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = t_end):
-            pass
-
-        # Check that momentum is as it should be in the interior
-
-        I = [[0, width/2.],
-             [length/2., width/2.],
-             [length, width/2.]]
-        
-        f = file_function(swwfile,
-                          quantities=['stage', 'xmomentum', 'ymomentum'],
-                          interpolation_points=I,
-                          verbose=False)
-        for t in range(t_end+1):
-            for i in range(3):
-                assert num.allclose(f(t, i), [1, 2, 0], atol=1.0e-6)
-            
-
-        # Check flows through the middle
-        for i in range(5):
-            x = length/2. + i*0.23674563 # Arbitrary
-            cross_section = [[x, 0], [x, width]]
-            time, Q = get_flow_through_cross_section(swwfile,
-                                                     cross_section,
-                                                     verbose=False)
-
-            assert num.allclose(Q, uh*width)
-
-
-       
-        # Try the same with partial lines
-        x = length/2.
-        for i in range(5):
-            start_point = [length/2., i*width/5.]
-            #print start_point
-                            
-            cross_section = [start_point, [length/2., width]]
-            time, Q = get_flow_through_cross_section(swwfile,
-                                                     cross_section,
-                                                     verbose=False)
-
-            #print i, Q, (width-start_point[1])
-            assert num.allclose(Q, uh*(width-start_point[1]))
-
-
-        # Verify no flow when line is parallel to flow
-        cross_section = [[length/2.-10, width/2.], [length/2.+10, width/2.]]
-        time, Q = get_flow_through_cross_section(swwfile,
-                                                 cross_section,
-                                                 verbose=False)
-
-        #print i, Q
-        assert num.allclose(Q, 0, atol=1.0e-5)        
-
-
-        # Try with lines on an angle (all flow still runs through here)
-        cross_section = [[length/2., 0], [length/2.+width, width]]
-        time, Q = get_flow_through_cross_section(swwfile,
-                                                 cross_section,
-                                                 verbose=False)
-
-        assert num.allclose(Q, uh*width)        
-        
-
-
-                                      
-    def test_get_flow_through_cross_section_with_geo(self):
-        """test_get_flow_through_cross_section(self):
-
-        Test that the total flow through a cross section can be
-        correctly obtained from an sww file.
-        
-        This test creates a flat bed with a known flow through it and tests
-        that the function correctly returns the expected flow.
-
-        The specifics are
-        u = 2 m/s
-        h = 2 m
-        w = 3 m (width of channel)
-
-        q = u*h*w = 12 m^3/s
-
-
-        This run tries it with georeferencing and with elevation = -1
-        
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        # Setup
-        from mesh_factory import rectangular
-
-        # Create basic mesh (20m x 3m)
-        width = 3
-        length = 20
-        t_end = 1
-        points, vertices, boundary = rectangular(length, width,
-                                                 length, width)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary,
-                        geo_reference = Geo_reference(56,308500,6189000))
-
-        domain.default_order = 2
-        domain.set_minimum_storable_height(0.01)
-
-        domain.set_name('flowtest')
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-
-        e = -1.0
-        w = 1.0
-        h = w-e
-        u = 2.0
-        uh = u*h
-
-        Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
-
-
-
-        
-        domain.set_quantity('elevation', e)
-        domain.set_quantity('stage', w)
-        domain.set_quantity('xmomentum', uh)
-        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = t_end):
-            pass
-
-        # Check that momentum is as it should be in the interior
-
-        I = [[0, width/2.],
-             [length/2., width/2.],
-             [length, width/2.]]
-        
-        I = domain.geo_reference.get_absolute(I)
-        f = file_function(swwfile,
-                          quantities=['stage', 'xmomentum', 'ymomentum'],
-                          interpolation_points=I,
-                          verbose=False)
-
-        for t in range(t_end+1):
-            for i in range(3):
-                #print i, t, f(t, i)            
-                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
-            
-
-        # Check flows through the middle
-        for i in range(5):
-            x = length/2. + i*0.23674563 # Arbitrary
-            cross_section = [[x, 0], [x, width]]
-
-            cross_section = domain.geo_reference.get_absolute(cross_section)            
-            time, Q = get_flow_through_cross_section(swwfile,
-                                                     cross_section,
-                                                     verbose=False)
-
-            assert num.allclose(Q, uh*width)
-
-
-            
-    def test_get_energy_through_cross_section(self):
-        """test_get_energy_through_cross_section(self):
-
-        Test that the specific and total energy through a cross section can be
-        correctly obtained from an sww file.
-        
-        This test creates a flat bed with a known flow through it and tests
-        that the function correctly returns the expected energies.
-
-        The specifics are
-        u = 2 m/s
-        h = 1 m
-        w = 3 m (width of channel)
-
-        q = u*h*w = 6 m^3/s
-        Es = h + 0.5*v*v/g  # Specific energy head [m]
-        Et = w + 0.5*v*v/g  # Total energy head [m]        
-
-
-        This test uses georeferencing
-        
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        # Setup
-        from mesh_factory import rectangular
-
-        # Create basic mesh (20m x 3m)
-        width = 3
-        length = 20
-        t_end = 1
-        points, vertices, boundary = rectangular(length, width,
-                                                 length, width)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary,
-                        geo_reference = Geo_reference(56,308500,6189000))
-
-        domain.default_order = 2
-        domain.set_minimum_storable_height(0.01)
-
-        domain.set_name('flowtest')
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-
-        e = -1.0
-        w = 1.0
-        h = w-e
-        u = 2.0
-        uh = u*h
-
-        Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
-
-        
-        domain.set_quantity('elevation', e)
-        domain.set_quantity('stage', w)
-        domain.set_quantity('xmomentum', uh)
-        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = t_end):
-            pass
-
-        # Check that momentum is as it should be in the interior
-
-        I = [[0, width/2.],
-             [length/2., width/2.],
-             [length, width/2.]]
-        
-        I = domain.geo_reference.get_absolute(I)
-        f = file_function(swwfile,
-                          quantities=['stage', 'xmomentum', 'ymomentum'],
-                          interpolation_points=I,
-                          verbose=False)
-
-        for t in range(t_end+1):
-            for i in range(3):
-                #print i, t, f(t, i)
-                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
-            
-
-        # Check energies through the middle
-        for i in range(5):
-            x = length/2. + i*0.23674563 # Arbitrary
-            cross_section = [[x, 0], [x, width]]
-
-            cross_section = domain.geo_reference.get_absolute(cross_section)            
-            
-            time, Es = get_energy_through_cross_section(swwfile,
-                                                       cross_section,
-                                                       kind='specific',
-                                                       verbose=False)
-            assert num.allclose(Es, h + 0.5*u*u/g)
-            
-            time, Et = get_energy_through_cross_section(swwfile,
-                                                        cross_section,
-                                                        kind='total',
-                                                        verbose=False)
-            assert num.allclose(Et, w + 0.5*u*u/g)            
-
-            
         
     def test_get_all_swwfiles(self):

trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -29 +29 @@
 from anuga.shallow_water.forcing import Rainfall, Wind_stress
 from anuga.shallow_water.forcing import Inflow, Cross_section
-from anuga.shallow_water.data_manager import get_flow_through_cross_section
+from anuga.shallow_water.sww_interrogate import get_flow_through_cross_section
 
 # boundary functions
@@ -1257 +1257 @@
         assert num.alltrue(wet_elevations < final_runup_height)
         assert num.allclose(wet_elevations, z)
-
-    def test_get_maximum_inundation_from_sww(self):
-        """test_get_maximum_inundation_from_sww(self)
-
-        Test of get_maximum_inundation_elevation()
-        and get_maximum_inundation_location() from data_manager.py
-
-        This is based on test_get_maximum_inundation_3(self) but works with the
-        stored results instead of with the internal data structure.
-
-        This test uses the underlying get_maximum_inundation_data for tests
-        """
-
-        from anuga.abstract_2d_finite_volumes.mesh_factory \
-                import rectangular_cross
-        from data_manager import get_maximum_inundation_elevation
-        from data_manager import get_maximum_inundation_location
-        from data_manager import get_maximum_inundation_data
-
-        initial_runup_height = -0.4
-        final_runup_height = -0.3
-
-        #--------------------------------------------------------------
-        # Setup computational domain
-        #--------------------------------------------------------------
-        N = 10
-        points, vertices, boundary = rectangular_cross(N, N)
-        domain = Domain(points, vertices, boundary)
-        domain.set_name('runup_test')
-        domain.set_maximum_allowed_speed(1.0)
-
-        # FIXME: This works better with old limiters so far
-        domain.tight_slope_limiters = 0
-
-        #--------------------------------------------------------------
-        # Setup initial conditions
-        #--------------------------------------------------------------
-        def topography(x, y):
-            return -x/2                             # linear bed slope
-
-        # Use function for elevation
-        domain.set_quantity('elevation', topography)
-        domain.set_quantity('friction', 0.)                # Zero friction
-        # Constant negative initial stage
-        domain.set_quantity('stage', initial_runup_height)
-
-        #--------------------------------------------------------------
-        # Setup boundary conditions
-        #--------------------------------------------------------------
-        Br = Reflective_boundary(domain)                       # Reflective wall
-        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
-
-        # All reflective to begin with (still water)
-        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
-
-        #--------------------------------------------------------------
-        # Test initial inundation height
-        #--------------------------------------------------------------
-        indices = domain.get_wet_elements()
-        z = domain.get_quantity('elevation').\
-                get_values(location='centroids', indices=indices)
-        assert num.alltrue(z < initial_runup_height)
-
-        q_ref = domain.get_maximum_inundation_elevation()
-        # First order accuracy
-        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
-
-        #--------------------------------------------------------------
-        # Let triangles adjust
-        #--------------------------------------------------------------
-        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
-            pass
-
-        #--------------------------------------------------------------
-        # Test inundation height again
-        #--------------------------------------------------------------
-        q_ref = domain.get_maximum_inundation_elevation()
-        q = get_maximum_inundation_elevation('runup_test.sww')
-        msg = 'We got %f, should have been %f' % (q, q_ref)
-        assert num.allclose(q, q_ref, rtol=1.0/N), msg
-
-        q = get_maximum_inundation_elevation('runup_test.sww')
-        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
-        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
-
-        # Test error condition if time interval is out
-        try:
-            q = get_maximum_inundation_elevation('runup_test.sww',
-                                                 time_interval=[2.0, 3.0])
-        except ValueError:
-            pass
-        else:
-            msg = 'should have caught wrong time interval'
-            raise Exception, msg
-
-        # Check correct time interval
-        q, loc = get_maximum_inundation_data('runup_test.sww',
-                                             time_interval=[0.0, 3.0])
-        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
-        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
-        assert num.allclose(-loc[0]/2, q)    # From topography formula
-
-        #--------------------------------------------------------------
-        # Update boundary to allow inflow
-        #--------------------------------------------------------------
-        domain.set_boundary({'right': Bd})
-
-        #--------------------------------------------------------------
-        # Evolve system through time
-        #--------------------------------------------------------------
-        q_max = None
-        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
-                               skip_initial_step = True):
-            q = domain.get_maximum_inundation_elevation()
-            if q > q_max:
-                q_max = q
-
-        #--------------------------------------------------------------
-        # Test inundation height again
-        #--------------------------------------------------------------
-        indices = domain.get_wet_elements()
-        z = domain.get_quantity('elevation').\
-                get_values(location='centroids', indices=indices)
-
-        assert num.alltrue(z < final_runup_height)
-
-        q = domain.get_maximum_inundation_elevation()
-        # First order accuracy
-        assert num.allclose(q, final_runup_height, rtol=1.0/N)
-
-        q, loc = get_maximum_inundation_data('runup_test.sww',
-                                             time_interval=[3.0, 3.0])
-        msg = 'We got %f, should have been %f' % (q, final_runup_height)
-        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
-        assert num.allclose(-loc[0]/2, q)    # From topography formula
-
-        q = get_maximum_inundation_elevation('runup_test.sww')
-        loc = get_maximum_inundation_location('runup_test.sww')
-        msg = 'We got %f, should have been %f' % (q, q_max)
-        assert num.allclose(q, q_max, rtol=1.0/N), msg
-        assert num.allclose(-loc[0]/2, q)    # From topography formula
-
-        q = get_maximum_inundation_elevation('runup_test.sww',
-                                             time_interval=[0, 3])
-        msg = 'We got %f, should have been %f' % (q, q_max)
-        assert num.allclose(q, q_max, rtol=1.0/N), msg
-
-        # Check polygon mode
-        # Runup region
-        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
-        q = get_maximum_inundation_elevation('runup_test.sww',
-                                             polygon = polygon,
-                                             time_interval=[0, 3])
-        msg = 'We got %f, should have been %f' % (q, q_max)
-        assert num.allclose(q, q_max, rtol=1.0/N), msg
-
-        # Offshore region
-        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
-        q, loc = get_maximum_inundation_data('runup_test.sww',
-                                             polygon = polygon,
-                                             time_interval=[0, 3])
-        msg = 'We got %f, should have been %f' % (q, -0.475)
-        assert num.allclose(q, -0.475, rtol=1.0/N), msg
-        assert is_inside_polygon(loc, polygon)
-        assert num.allclose(-loc[0]/2, q)    # From topography formula
-
-        # Dry region
-        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]]
-        q, loc = get_maximum_inundation_data('runup_test.sww',
-                                             polygon = polygon,
-                                             time_interval=[0, 3])
-        msg = 'We got %s, should have been None' % (q)
-        assert q is None, msg
-        msg = 'We got %s, should have been None' % (loc)
-        assert loc is None, msg
-
-        # Check what happens if no time point is within interval
-        try:
-            q = get_maximum_inundation_elevation('runup_test.sww',
-                                                 time_interval=[2.75, 2.75])
-        except AssertionError:
-            pass
-        else:
-            msg = 'Time interval should have raised an exception'
-            raise Exception, msg
-
-        # Cleanup
-        try:
-            os.remove(domain.get_name() + '.sww')
-        except:
-            pass
-            #FIXME(Ole): Windows won't allow removal of this
 
     def test_get_flow_through_cross_section_with_geo(self):
@@ -6085 +5893 @@
         from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
         from anuga.shallow_water.shallow_water_domain import Rainfall
-        from anuga.shallow_water.data_manager import get_flow_through_cross_section
+        from anuga.shallow_water.sww_interrogate import get_flow_through_cross_section
 
         #----------------------------------------------------------------------