Changeset 6689 for branches/numpy/anuga

Timestamp:

Apr 1, 2009, 3:19:07 PM (16 years ago)

Author:

rwilson

Message:

Back-merge from Numeric trunk to numpy branch.

Location:

branches/numpy/anuga

Files:

• abstract_2d_finite_volumes/general_mesh.py (modified) (1 diff)
• abstract_2d_finite_volumes/pmesh2domain.py (modified) (2 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (1 diff)
• abstract_2d_finite_volumes/test_general_mesh.py (modified) (2 diffs)
• abstract_2d_finite_volumes/test_region.py (modified) (1 diff)
• abstract_2d_finite_volumes/util.py (modified) (1 diff)
• advection/test_advection.py (modified) (1 diff)
• caching/test_caching.py (modified) (1 diff)
• coordinate_transforms/geo_reference.py (modified) (1 diff)
• culvert_flows/culvert_routines.py (modified) (5 diffs)
• culvert_flows/test_culvert_routines.py (modified) (3 diffs)
• fit_interpolate/fit.py (modified) (5 diffs)
• fit_interpolate/interpolate.py (modified) (2 diffs)
• fit_interpolate/test_interpolate.py (modified) (3 diffs)
• fit_interpolate/test_search_functions.py (modified) (1 diff)
• geospatial_data/geospatial_data.py (modified) (2 diffs)
• pmesh/mesh_interface.py (modified) (1 diff)
• shallow_water/data_manager.py (modified) (17 diffs)
• shallow_water/shallow_water_domain.py (modified) (7 diffs)
• shallow_water/test_data_manager.py (modified) (20 diffs)
• shallow_water/test_points_large.csv (added)
• shallow_water/test_points_large.lic (added)
• shallow_water/test_points_small.csv (added)
• shallow_water/test_points_small.lic (added)
• shallow_water/test_shallow_water_domain.py (modified) (16 diffs)
• utilities/log.py (modified) (8 diffs)
• utilities/polygon.py (modified) (2 diffs)
• utilities/system_tools.py (modified) (7 diffs)
• utilities/test_polygon.py (modified) (1 diff)
• utilities/test_system_tools.py (modified) (3 diffs)
• utilities/util_ext.h (modified) (1 diff)

branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -221 +221 @@
 
         return self.normals[i, 2*j:2*j+2]
+        
+    def get_edgelength(self, i, j):
+        """Return length of j'th edge of the i'th triangle.
+        Return value is the numeric array slice [x, y]
+        """
+        return self.edgelengths[i, j]
+                
 
     def get_number_of_nodes(self):

branches/numpy/anuga/abstract_2d_finite_volumes/pmesh2domain.py

@@ -165 +165 @@
                      verbose=False):
     """Convert a pmesh file or a pmesh mesh instance to a bunch of lists
-    that can be used to instanciate a domain object.
+    that can be used to instantiate a domain object.
     """
 
@@ -180 +180 @@
     volumes = mesh_dict['triangles']
     vertex_quantity_dict = {}
+
+    # num.transpose(None) gives scalar array of value None
     point_atts = mesh_dict['vertex_attributes']
-    point_titles  = mesh_dict['vertex_attribute_titles']
-    geo_reference  = mesh_dict['geo_reference']
+
+    point_titles = mesh_dict['vertex_attribute_titles']
+    geo_reference = mesh_dict['geo_reference']
     if point_atts is not None:
-        point_atts = num.transpose(mesh_dict['vertex_attributes'])
+        point_atts = num.transpose(point_atts)
         for quantity, value_vector in map(None, point_titles, point_atts):
             vertex_quantity_dict[quantity] = value_vector

branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

@@ -1188 +1188 @@
         msg = 'Indices must be a list or None'
         assert indices is None or isinstance(indices, (list, num.ndarray)), msg
-#        assert type(indices) in [types.ListType, types.NoneType,
-#                                 num.ndarray], \
-#                   'Indices must be a list or None'     #??#
 
         if location == 'centroids':

branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -132 +132 @@
         domain = General_mesh(nodes, triangles)
 
-#        value = [7]
         msg = ("domain.get_triangles()=\n%s\nshould be the same as "
                "'triangles'=\n%s"
@@ -354 +353 @@
                               nodes, triangles, geo_reference=geo)
 
-    def test_raw_change_points_geo_ref(self):
-        x0 = 1000.0
-        y0 = 2000.0
-        geo = Geo_reference(56, x0, y0)
-        
-        
-
-
-#-------------------------------------------------------------
+################################################################################
 
 if __name__ == "__main__":

branches/numpy/anuga/abstract_2d_finite_volumes/test_region.py

@@ -68 +68 @@
                     [ 0.07,  0.07,  0.07],
                     [ 0.07,  0.07,  0.07],
-                    [ 1.0,  1.0,  1.0],
-                    [ 1.0,  1.0,  1.0]]
+                    [ 1.0,   1.0,   1.0],
+                    [ 1.0,   1.0,   1.0]]
         msg = ("\ndomain.quantities['friction']=%s\nexpected value=%s"
                % (str(domain.quantities['friction'].get_values()),

branches/numpy/anuga/abstract_2d_finite_volumes/util.py

@@ -2122 +2122 @@
                 #add tide to stage if provided
                 if quantity == 'stage':
-                     quantity_value[quantity] = num.array(quantity_value[quantity], num.float) \
-                                                          + directory_add_tide
+                     quantity_value[quantity] = num.array(quantity_value[quantity],
+                                                          num.float) + directory_add_tide
 
                 #condition to find max and mins for all the plots

branches/numpy/anuga/advection/test_advection.py

@@ -142 +142 @@
 
         X = domain.quantities['stage'].explicit_update
-        #        print 'X=%s' % str(X)
         assert X[0] == -X[1]
 

branches/numpy/anuga/caching/test_caching.py

@@ -390 +390 @@
         if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
           # 64 bit hash values
-          f1hash = 1914027059797211698
-          f2hash = 1914027059807087171
+          f1hash = 7079146893884768701
+          f2hash = -6995306676314913340
+
+          # Prior to cluster upgrades Feb 2009
+          #f1hash = 1914027059797211698 
+          #f2hash = 1914027059807087171
         else:
+          # 32 bit hash values
           f1hash = -758136387
           f2hash = -11221564     

branches/numpy/anuga/coordinate_transforms/geo_reference.py

@@ -258 +258 @@
         is_list = isinstance(points, list)
 
-        #points = ensure_numeric(copy.copy(points), num.float)
         points = ensure_numeric(points, num.float)
 

branches/numpy/anuga/culvert_flows/culvert_routines.py

@@ -28 +28 @@
                                    manning=0.0,
                                    sum_loss=0.0,
+                                   max_velocity=10.0,
                                    log_filename=None):
 
@@ -47 +48 @@
 
 
-    if inlet_depth > 0.01:
+    if inlet_depth > 0.1: #this value was 0.01:
         # Water has risen above inlet
         
@@ -116 +117 @@
 
                 # Outlet control velocity using tail water
-                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
+                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
                 Q_outlet_tailwater = flow_area * culvert_velocity
                 
@@ -183 +184 @@
 
                 # Outlet control velocity using tail water
-                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
+                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
                 Q_outlet_tailwater = flow_area * culvert_velocity
 
@@ -216 +217 @@
     else: # inlet_depth < 0.01:
         Q = barrel_velocity = outlet_culvert_depth = 0.0
+    # Temporary flow limit
+    if barrel_velocity > max_velocity:
+        if log_filename is not None:                            
+            s = 'Barrel velocity was reduced from = %f m/s to %f m/s' % (barrel_velocity, max_velocity)
+            log_to_file(log_filename, s)
+        
+        barrel_velocity = max_velocity
+        Q = flow_area * barrel_velocity
+        
+        
+
         
     return Q, barrel_velocity, outlet_culvert_depth 

branches/numpy/anuga/culvert_flows/test_culvert_routines.py

@@ -19 +19 @@
 
 
-    def NOtest_boyd_1(self):
+    def test_boyd_1(self):
         """test_boyd_1
         
@@ -29 +29 @@
         culvert_slope=0.1  # Downward
 
-        inlet_depth=0.1
-        outlet_depth=0.09
+        inlet_depth=2.0
+        outlet_depth=0.0
 
         culvert_length=4.0
@@ -59 +59 @@
                                                  sum_loss)
         
-        print Q, v, d
-        assert num.allclose(Q, 0.1)
-        assert num.allclose(v, 0.93)
-        assert num.allclose(d, 0.09)
+        #print Q, v, d
+        assert num.allclose(Q, 3.118, rtol=1.0e-3)
         
 
+        #assert num.allclose(v, 0.93)
+        #assert num.allclose(d, 0.0)
+        
+
+    def test_boyd_2(self):
+        """test_boyd_2
+        
+        This tests the Boyd routine with data obtained from ??? by Petar Milevski    
+        """
+        # FIXME(Ole): This test fails (20 Feb 2009)
+
+        g=9.81
+        culvert_slope=0.1  # Downward
+
+        inlet_depth=0.2
+        outlet_depth=0.0
+
+        culvert_length=4.0
+        culvert_width=1.2
+        culvert_height=0.75
+
+        culvert_type='box'
+        manning=0.013
+        sum_loss=0.0
+
+        inlet_specific_energy=inlet_depth #+0.5*v**2/g 
+        z_in = 0.0
+        z_out = -culvert_length*culvert_slope/100
+        E_in = z_in+inlet_depth # + 
+        E_out = z_out+outlet_depth # +
+        delta_total_energy = E_in-E_out
+
+        Q, v, d = boyd_generalised_culvert_model(inlet_depth, 
+                                                 outlet_depth,
+                                                 inlet_specific_energy, 
+                                                 delta_total_energy, 
+                                                 g,
+                                                 culvert_length,
+                                                 culvert_width,
+                                                 culvert_height,
+                                                 culvert_type,
+                                                 manning,
+                                                 sum_loss)
+        
+        #print Q, v, d
+        #assert num.allclose(Q, 0.185, rtol=1.0e-3)
+        #assert num.allclose(v, 0.93)
+        #assert num.allclose(d, 0.0)
+        
     
                

branches/numpy/anuga/fit_interpolate/fit.py

@@ -260 +260 @@
             # The memory damage has been done by now.
         else:
-             AtA = self.AtA #Did this for speed, did ~nothing
+             AtA = self.AtA # Did this for speed, did ~nothing
         self.point_count += point_coordinates.shape[0]
 
@@ -267 +267 @@
                                         self.mesh_boundary_polygon,
                                         closed=True,
-                                        verbose=False) # Too much output if True
-
+                                        verbose=False) # Suppress output
         
         n = len(inside_indices)
@@ -283 +282 @@
             
             if element_found is True:
-                j0 = triangles[k,0] #Global vertex id for sigma0
-                j1 = triangles[k,1] #Global vertex id for sigma1
-                j2 = triangles[k,2] #Global vertex id for sigma2
+                j0 = triangles[k,0] # Global vertex id for sigma0
+                j1 = triangles[k,1] # Global vertex id for sigma1
+                j2 = triangles[k,2] # Global vertex id for sigma2
 
                 sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
@@ -304 +303 @@
                                          self.mesh_boundary_polygon,
                                          closed=True,
-                                         verbose=False) # Too much output if True                
+                                         verbose=False) # Suppress output
                 msg = 'Point (%f, %f) is not inside mesh boundary' % tuple(x)
                 assert flag is True, msg                
@@ -318 +317 @@
                 raise RuntimeError, msg
 
-            self.AtA = AtA
+                
+        self.AtA = AtA
 
         

branches/numpy/anuga/fit_interpolate/interpolate.py

@@ -835 +835 @@
                 raise msg
 
+            # Ensure 'mesh_boundary_polygon' is defined
+            mesh_boundary_polygon = None
+            
             if triangles is not None and vertex_coordinates is not None:
                 # Check that all interpolation points fall within
@@ -894 +897 @@
                                 'requires the latter.')
 
-            # Plot boundary and interpolation points
-            if verbose is True:
+            # Plot boundary and interpolation points,
+            # but only if if 'mesh_boundary_polygon' has data.
+            if verbose is True and mesh_boundary_polygon is not None:
                 import sys
                 if sys.platform == 'win32':

branches/numpy/anuga/fit_interpolate/test_interpolate.py

@@ -919 +919 @@
         triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
 
-        point_coords_absolute = [[-2.0, 2.0],
-                                 [-1.0, 1.0],
-                                 [0.0, 2.0],
-                                 [1.0, 1.0],
-                                 [2.0, 1.0],
-                                 [0.0, 0.0],
-                                 [1.0, 0.0],
-                                 [0.0, -1.0],
+        point_coords_absolute = [[-2.0,  2.0],
+                                 [-1.0,  1.0],
+                                 [ 0.0,  2.0],
+                                 [ 1.0,  1.0],
+                                 [ 2.0,  1.0],
+                                 [ 0.0,  0.0],
+                                 [ 1.0,  0.0],
+                                 [ 0.0, -1.0],
                                  [-0.2, -0.5],
                                  [-0.9, -1.5],
-                                 [0.5, -1.9],
-                                 [3.0, 1.0]]
+                                 [ 0.5, -1.9],
+                                 [ 3.0,  1.0]]
 
         geo = Geo_reference(57, 100, 500)
@@ -967 +967 @@
 
 
-        point_coords = [[-2.0, 2.0],
-                        [-1.0, 1.0],
-                        [0.0, 2.0],
-                        [1.0, 1.0],
-                        [2.0, 1.0],
-                        [0.0, 0.0],
-                        [1.0, 0.0],
-                        [0.0, -1.0],
+        point_coords = [[-2.0,  2.0],
+                        [-1.0,  1.0],
+                        [ 0.0,  2.0],
+                        [ 1.0,  1.0],
+                        [ 2.0,  1.0],
+                        [ 0.0,  0.0],
+                        [ 1.0,  0.0],
+                        [ 0.0, -1.0],
                         [-0.2, -0.5],
                         [-0.9, -1.5],
-                        [0.5, -1.9],
-                        [3.0, 1.0]]
+                        [ 0.5, -1.9],
+                        [ 3.0,  1.0]]
 
         interp = Interpolate(vertices, triangles)
@@ -1030 +1030 @@
 
 
-        point_coords = [[-2.0, 2.0],
-                        [-1.0, 1.0],
-                        [0.0, 2.0],
-                        [1.0, 1.0],
-                        [2.0, 1.0],
-                        [0.0, 0.0],
-                        [1.0, 0.0],
-                        [0.0, -1.0],
+        point_coords = [[-2.0,  2.0],
+                        [-1.0,  1.0],
+                        [ 0.0,  2.0],
+                        [ 1.0,  1.0],
+                        [ 2.0,  1.0],
+                        [ 0.0,  0.0],
+                        [ 1.0,  0.0],
+                        [ 0.0, -1.0],
                         [-0.2, -0.5],
                         [-0.9, -1.5],
-                        [0.5, -1.9],
-                        [3.0, 1.0]]
+                        [ 0.5, -1.9],
+                        [ 3.0,  1.0]]
 
         interp = Interpolate(vertices, triangles)

branches/numpy/anuga/fit_interpolate/test_search_functions.py

@@ -113 +113 @@
                 if k >= 0:
                     V = mesh.get_vertex_coordinates(k) # nodes for triangle k
+                    assert is_inside_triangle(x, V, closed=True)
                     assert is_inside_polygon(x, V)
                     assert found is True

branches/numpy/anuga/geospatial_data/geospatial_data.py

@@ -11 +11 @@
 from RandomArray import randint, seed, get_seed
 from copy import deepcopy
+import copy
 
 from Scientific.IO.NetCDF import NetCDFFile
@@ -1418 +1419 @@
         assert geo_reference == None, msg
     else:
-        points = ensure_numeric(points, num.float)
+        points = ensure_numeric(copy.copy(points), num.float)
 
     # Sort of geo_reference and convert points

branches/numpy/anuga/pmesh/mesh_interface.py

@@ -147 +147 @@
                     found = True
             if found is False:
-                msg = 'Segment %d was not asigned a boundary_tag' %i
-                raise Exception, msg
+                msg = 'Segment %d was not assigned a boundary_tag.' % i
+                msg +=  'Default tag "exterior" will be assigned to missing segment'
+                #raise Exception, msg
+                # Fixme: Use proper Python warning
+                if verbose: print 'WARNING: ', msg
                 
 

branches/numpy/anuga/shallow_water/data_manager.py

@@ -91 +91 @@
 from anuga.utilities.polygon import intersection
 
+from anuga.utilities.system_tools import get_vars_in_expression
+
+
+# Default block size for sww2dem()
+DEFAULT_BLOCK_SIZE = 10000
 
 ######
@@ -2106 +2111 @@
             origin=None,
             datum='WGS84',
-            format='ers'):
+            format='ers',
+            block_size=None):
     """Read SWW file and convert to Digitial Elevation model format
     (.asc or .ers)
@@ -2149 +2155 @@
 
     format can be either 'asc' or 'ers'
+    block_size - sets the number of slices along the non-time axis to
+                 process in one block.
     """
 
@@ -2179 +2187 @@
         number_of_decimal_places = 3
 
+    if block_size is None:
+        block_size = DEFAULT_BLOCK_SIZE
+
+    # Read SWW file
     swwfile = basename_in + '.sww'
     demfile = basename_out + '.' + format
@@ -2255 +2267 @@
             if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
 
-    # Get quantity and reduce if applicable
-    if verbose: print 'Processing quantity %s' %quantity
-
-    # Turn NetCDF objects into numeric arrays
-    try:
-        q = fid.variables[quantity][:]
-    except:
-        quantity_dict = {}
-        for name in fid.variables.keys():
-            quantity_dict[name] = fid.variables[name][:]
-        #Convert quantity expression to quantities found in sww file
-        q = apply_expression_to_dictionary(quantity, quantity_dict)
-
-    if len(q.shape) == 2:
-        #q has a time component, must be reduced alongthe temporal dimension
-        if verbose: print 'Reducing quantity %s' %quantity
-        q_reduced = num.zeros(number_of_points, num.float)
-
-        if timestep is not None:
-            for k in range(number_of_points):
-                q_reduced[k] = q[timestep,k]
-        else:
-            for k in range(number_of_points):
-                q_reduced[k] = reduction( q[:,k] )
-
-        q = q_reduced
-
+    # Get the variables in the supplied expression.
+    # This may throw a SyntaxError exception.
+    var_list = get_vars_in_expression(quantity)
+
+    # Check that we have the required variables in the SWW file.
+    missing_vars = []
+    for name in var_list:
+        try:
+            _ = fid.variables[name]
+        except:
+            missing_vars.append(name)
+    if missing_vars:
+        msg = ("In expression '%s', variables %s are not in the SWW file '%s'"
+               % (quantity, swwfile))
+        raise Exception, msg
+
+    # Create result array and start filling, block by block.
+    result = num.zeros(number_of_points, num.float)
+
+    for start_slice in xrange(0, number_of_points, block_size):
+        # limit slice size to array end if at last block
+        end_slice = min(start_slice + block_size, number_of_points)
+        
+        # get slices of all required variables
+        q_dict = {}
+        for name in var_list:
+            # check if variable has time axis
+            if len(fid.variables[name].shape) == 2:
+                q_dict[name] = fid.variables[name][:,start_slice:end_slice]
+            else:       # no time axis
+                q_dict[name] = fid.variables[name][start_slice:end_slice]
+
+        # Evaluate expression with quantities found in SWW file
+        res = apply_expression_to_dictionary(quantity, q_dict)
+
+        if len(res.shape) == 2:
+            new_res = num.zeros(res.shape[1], num.float)
+            for k in xrange(res.shape[1]):
+                new_res[k] = reduction(res[:,k])
+            res = new_res
+
+        result[start_slice:end_slice] = res
+                                    
     #Post condition: Now q has dimension: number_of_points
-    assert len(q.shape) == 1
-    assert q.shape[0] == number_of_points
+    assert len(result.shape) == 1
+    assert result.shape[0] == number_of_points
 
     if verbose:
         print 'Processed values for %s are in [%f, %f]' % \
-              (quantity, min(q), max(q))
+              (quantity, min(result), max(result))
 
     #Create grid and update xll/yll corner and x,y
@@ -2316 +2344 @@
     assert xmax >= xmin, msg
 
-    msg = 'yax must be greater than or equal to xmin.\n'
+    msg = 'ymax must be greater than or equal to xmin.\n'
     msg += 'I got ymin = %f, ymax = %f' %(ymin, ymax)
     assert ymax >= ymin, msg
@@ -2361 +2389 @@
     #Interpolate using quantity values
     if verbose: print 'Interpolating'
-    grid_values = interp.interpolate(q, grid_points).flatten()
+    grid_values = interp.interpolate(result, grid_points).flatten()
 
     if verbose:
@@ -2672 +2700 @@
 
     if verbose:
-        print 'Interpolated values are in [%f, %f]' \
-              % (num.min(interpolated_values), num.max(interpolated_values))
+        print ('Interpolated values are in [%f, %f]'
+               % (num.min(interpolated_values), num.max(interpolated_values)))
 
     # Assign NODATA_value to all points outside bounding polygon
@@ -5475 +5503 @@
             origin=None,
             zone=None,
+            central_meridian=None,            
             mean_stage=0.0,
             zscale=1.0,
@@ -5514 +5543 @@
     output:
       basename_out: name of sts file in which mux2 data is stored.
+      
+      
+      
+    NOTE: South is positive in mux files so sign of y-component of velocity is reverted
     """
 
@@ -5687 +5720 @@
     # Check zone boundaries
     if zone is None:
-        refzone, _, _ = redfearn(latitudes[0], longitudes[0])
+        refzone, _, _ = redfearn(latitudes[0], longitudes[0],
+                                 central_meridian=central_meridian)
     else:
         refzone = zone
@@ -5694 +5728 @@
 
     for i in range(number_of_points):
-        zone, easting, northing = redfearn(latitudes[i], longitudes[i],
-                                           zone=zone)
+        computed_zone, easting, northing = redfearn(latitudes[i], longitudes[i],
+                                                    zone=zone,
+                                                    central_meridian=central_meridian)
         x[i] = easting
         y[i] = northing
-        if zone != refzone:
+        if computed_zone != refzone:
             msg = 'All sts gauges need to be in the same zone. \n'
             msg += 'offending gauge:Zone %d,%.4f, %4f\n' \
-                   % (zone, easting, northing)
+                   % (computed_zone, easting, northing)
             msg += 'previous gauge:Zone %d,%.4f, %4f' \
                    % (old_zone, old_easting, old_northing)
             raise Exception, msg
-        old_zone = zone
+        old_zone = computed_zone
         old_easting = easting
         old_northing = northing
@@ -5744 +5779 @@
 
             xmomentum[j,i] = ua * h
-            ymomentum[j,i] = va * h
+            ymomentum[j,i] = -va * h # South is positive in mux files
+
 
     outfile.close()
@@ -5863 +5899 @@
         # This is being used to seperate one number from a list.
         # what it is actually doing is sorting lists from numeric arrays.
-        if type(times) is list or isinstance(times, num.ndarray):
+        if isinstance(times, (list, num.ndarray)):
             number_of_times = len(times)
             times = ensure_numeric(times)
@@ -5930 +5966 @@
             #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
 
-        if type(times) is list or isinstance(times, num.ndarray):
+        if isinstance(times, (list, num.ndarray)):
             outfile.variables['time'][:] = times    #Store time relative
 
@@ -6265 +6301 @@
         # This is being used to seperate one number from a list.
         # what it is actually doing is sorting lists from numeric arrays.
-        if type(times) is list or isinstance(times, num.ndarray):
+        if isinstance(times, (list, num.ndarray)):
             number_of_times = len(times)
             times = ensure_numeric(times)
@@ -6312 +6348 @@
             outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
 
-        if type(times) is list or isinstance(times, num.ndarray):
+        if isinstance(times, (list, num.ndarray)):
             outfile.variables['time'][:] = times    #Store time relative
 

branches/numpy/anuga/shallow_water/shallow_water_domain.py

@@ -812 +812 @@
 
         return msg
-
+       
+        
+
+    def compute_boundary_flows(self):
+        """Compute boundary flows at current timestep.
+                        
+        Quantities computed are:
+           Total inflow across boundary
+           Total outflow across boundary
+           Flow across each tagged boundary segment
+        """
+                
+        # Run through boundary array and compute for each segment
+        # the normal momentum ((uh, vh) dot normal) times segment length.
+        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
+                        
+        # Compute flows along boundary
+        
+        uh = self.get_quantity('xmomentum').get_values(location='edges')
+        vh = self.get_quantity('ymomentum').get_values(location='edges')        
+        
+        # Loop through edges that lie on the boundary and calculate 
+        # flows
+        boundary_flows = {}
+        total_boundary_inflow = 0.0
+        total_boundary_outflow = 0.0
+        for vol_id, edge_id in self.boundary:
+            # Compute normal flow across edge. Since normal vector points
+            # away from triangle, a positive sign means that water 
+            # flows *out* from this triangle.
+             
+            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
+            normal = self.mesh.get_normal(vol_id, edge_id)
+            length = self.mesh.get_edgelength(vol_id, edge_id)            
+            normal_flow = num.dot(momentum, normal)*length
+            
+            # Reverse sign so that + is taken to mean inflow
+            # and - means outflow. This is more intuitive.
+            edge_flow = -normal_flow
+            
+            # Tally up inflows and outflows separately
+            if edge_flow > 0:
+                # Flow is inflow      
+                total_boundary_inflow += edge_flow                                  
+            else:
+                # Flow is outflow
+                total_boundary_outflow += edge_flow    
+
+            # Tally up flows by boundary tag
+            tag = self.boundary[(vol_id, edge_id)]
+            
+            if tag not in boundary_flows:
+                boundary_flows[tag] = 0.0
+            boundary_flows[tag] += edge_flow
+            
+                
+        return boundary_flows, total_boundary_inflow, total_boundary_outflow
+        
+
+    def compute_forcing_flows(self):
+        """Compute flows in and out of domain due to forcing terms.
+                        
+        Quantities computed are:
+                
+        
+           Total inflow through forcing terms
+           Total outflow through forcing terms
+           Current total volume in domain        
+
+        """
+
+        #FIXME(Ole): We need to separate what part of explicit_update was 
+        # due to the normal flux calculations and what is due to forcing terms.
+        
+        pass
+                        
+        
+    def compute_total_volume(self):
+        """Compute total volume (m^3) of water in entire domain
+        """
+        
+        area = self.mesh.get_areas()
+        volume = 0.0
+        
+        stage = self.get_quantity('stage').get_values(location='centroids')
+        elevation = self.get_quantity('elevation').get_values(location='centroids')        
+        depth = stage-elevation
+        
+        #print 'z', elevation
+        #print 'w', stage
+        #print 'h', depth
+        return num.sum(depth*area)
+        
+        
+    def volumetric_balance_statistics(self):                
+        """Create volumetric balance report suitable for printing or logging.
+        """
+        
+        boundary_flows, total_boundary_inflow, total_boundary_outflow = self.compute_boundary_flows() 
+        
+        s = '---------------------------\n'        
+        s += 'Volumetric balance report:\n'
+        s += '--------------------------\n'
+        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
+        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow        
+        s += 'Net boundary flow by tags [m^3/s]\n'
+        for tag in boundary_flows:
+            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
+        
+        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 
+        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
+        
+        # The go through explicit forcing update and record the rate of change for stage and 
+        # record into forcing_inflow and forcing_outflow. Finally compute integral 
+        # of depth to obtain total volume of domain.
+        
+        # FIXME(Ole): This part is not yet done.                
+        
+        return s        
+           
 ################################################################################
 # End of class Shallow Water Domain
@@ -1386 +1505 @@
         domain: pointer to shallow water domain for which the boundary applies
         mean_stage: The mean water level which will be added to stage derived
-                    from the sww file
+                    from the boundary condition
         time_thinning: Will set how many time steps from the sww file read in
                        will be interpolated to the boundary. For example if
@@ -1402 +1521 @@
                           that available in the underlying data.
 
+                          Note that mean_stage will also be added to this.
+                                                
         use_cache:
         verbose:
@@ -1574 +1695 @@
                 xmom_update[k] += S*uh[k]
                 ymom_update[k] += S*vh[k]
+
+def depth_dependent_friction(domain, default_friction,
+                             surface_roughness_data,
+                             verbose=False):
+    """Returns an array of friction values for each wet element adjusted for depth.
+
+    Inputs:
+        domain - computational domain object
+        default_friction - depth independent bottom friction
+        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
+        friction region.
+
+    Outputs:
+        wet_friction - Array that can be used directly to update friction as follows:
+                       domain.set_quantity('friction', wet_friction)
+
+        
+        
+    """
+
+    import Numeric as num
+    
+    # Create a temp array to store updated depth dependent friction for wet elements
+    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
+    N=len(domain)
+    wet_friction    = num.zeros(N, num.Float)
+    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
+    
+    
+    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
+    # Recompute depth as vector  
+    d = depth.get_values(location='centroids')
+ 
+    # rebuild the 'friction' values adjusted for depth at this instant
+    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
+        
+        # Get roughness data for each element
+        n0 = float(surface_roughness_data[i,0])
+        d1 = float(surface_roughness_data[i,1])
+        n1 = float(surface_roughness_data[i,2])
+        d2 = float(surface_roughness_data[i,3])
+        n2 = float(surface_roughness_data[i,4])
+        
+        
+        # Recompute friction values from depth for this element 
+               
+        if d[i]   <= d1: 
+            depth_dependent_friction = n1
+        elif d[i] >= d2:
+            depth_dependent_friction = n2
+        else:
+            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
+            
+        # check sanity of result
+        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
+            print model_data.basename+' >>>> WARNING: computed depth_dependent friction out of range ddf,n1,n2 ', depth_dependent_friction, n1,n2
+        
+        # update depth dependent friction  for that wet element
+        wet_friction[i] = depth_dependent_friction
+        
+    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
+    
+    if verbose :
+        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
+        n_min=min(nvals)
+        n_max=max(nvals)
+        
+        print "         ++++ calculate_depth_dependent_friction  - Updated friction - range  %7.3f to %7.3f" %(n_min,n_max)
+    
+    return wet_friction
+
 
 ################################################################################
@@ -1901 +2093 @@
                 assert is_inside_polygon(point, bounding_polygon), msg
 
-            # Compute area and check that it is greater than 0
-            self.exchange_area = polygon_area(self.polygon)
-
-            msg = 'Polygon %s in forcing term' % str(self.polygon)
-            msg += ' has area = %f' % self.exchange_area
-            assert self.exchange_area > 0.0
-
         # Pointer to update vector
         self.update = domain.quantities[self.quantity_name].explicit_update
@@ -1931 +2116 @@
         if self.polygon is not None:
             # Inlet is polygon
-            inlet_region = ('polygon=\n%s, area=%f m^2' %
-                            (self.polygon, self.exchange_area))
-
+            inlet_region = 'polygon=%s' % (self.polygon) 
             self.exchange_indices = inside_polygon(points, self.polygon)
 
-        if self.exchange_indices is not None:
+        if self.exchange_indices is None:
+            self.exchange_area = polygon_area(bounding_polygon)
+        else:    
             if len(self.exchange_indices) == 0:
                 msg = 'No triangles have been identified in '
@@ -1942 +2127 @@
                 raise Exception, msg
 
+            # Compute exchange area as the sum of areas of triangles identified
+            # by circle or polygon
+            self.exchange_area = 0.0
+            for i in self.exchange_indices:
+                self.exchange_area += domain.areas[i]
+            
+
+        msg = 'Exchange area in forcing term'
+        msg += ' has area = %f' %self.exchange_area
+        assert self.exchange_area > 0.0            
+            
+                
+
+            
         # Check and store default_rate
         msg = ('Keyword argument default_rate must be either None '

branches/numpy/anuga/shallow_water/test_data_manager.py

@@ -1221 +1221 @@
 
     def test_sww2dem_asc_elevation_depth(self):
-        """
-        test_sww2dem_asc_elevation_depth(self):
+        """test_sww2dem_asc_elevation_depth
+        
         Test that sww information can be converted correctly to asc/prj
         format readable by e.g. ArcView
+        
+        Also check geo_reference is correct
         """
 
@@ -1249 +1251 @@
         sww.store_timestep()
 
-        #self.domain.tight_slope_limiters = 1
 
         self.domain.evolve_to_end(finaltime = 0.01)
@@ -2000 +2001 @@
                 number_of_decimal_places = 9,
                 verbose = self.verbose,
-                format = 'asc')
+                format = 'asc',
+                block_size=2)
 
 
@@ -5064 +5066 @@
         #print "longitudes_news",longitudes_news
 
-        self.failUnless(latitudes_new == [-35, -40, -45] and \
+        self.failUnless(latitudes_new == [-35, -40, -45] and
                         longitudes_news == [148, 149,150],
                          'failed')
@@ -5519 +5521 @@
         na and va quantities will be the Easting values.
         Depth and ua will be the Northing value.
+        
+        The mux file format has south as positive so
+        this function will swap the sign for va.  
         """
 
@@ -6036 +6041 @@
                 quantities_init[2].append(num.ones(time_step_count,num.float)*this_va) #
             else:
-                quantities_init[2].append(va[i])            
+                quantities_init[2].append(-va[i]) # South is negative in MUX
 
         file_handle, base_name = tempfile.mkstemp("")
@@ -6158 +6163 @@
         assert num.allclose(xvelocity,ua)
         msg='incorrect gauge va time series returned'
-        assert num.allclose(yvelocity,va)
+        assert num.allclose(yvelocity, -va)
 
     def test_urs2sts_read_mux2_pyII(self):
@@ -6216 +6221 @@
         assert num.allclose(xvelocity,ua)
         msg='incorrect gauge va time series returned'
-        assert num.allclose(yvelocity,va)
+        assert num.allclose(yvelocity, -va) # South is positive in MUX
 
     def test_urs2sts_read_mux2_pyIII(self):
@@ -6288 +6293 @@
         assert num.allclose(xvelocity,ua)
         msg='incorrect gauge va time series returned'
-        assert num.allclose(yvelocity,va)
+        assert num.allclose(yvelocity, -va) # South is positive in mux
         
 
@@ -6390 +6395 @@
                 if j == 0: assert num.allclose(data[i][:parameters_index], ha0[permutation[i], :])
                 if j == 1: assert num.allclose(data[i][:parameters_index], ua0[permutation[i], :])
-                if j == 2: assert num.allclose(data[i][:parameters_index], va0[permutation[i], :])
+                if j == 2: assert num.allclose(data[i][:parameters_index], -va0[permutation[i], :])
         
 
@@ -6596 +6601 @@
                         x = unpack('f',f.read(4))[0]
                         #print time, x, q_time[time, point_i]
+                        if q == 'VA': x = -x # South is positive in MUX
                         assert abs(q_time[time, point_i]-x)<epsilon
 
@@ -6669 +6675 @@
                     #print 'v ', data[i][:parameters_index][8]                    
                 
-                    assert num.allclose(data[i][:parameters_index], va1[permutation[i], :])
+                    # South is positive in MUX
+                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
                     
         
@@ -6787 +6794 @@
         os.remove(sts_file)
 
-    def test_urs2sts_nonstandard_projection(self):
+    def test_urs2sts_nonstandard_meridian(self):
         """
-        Test single source
+        Test single source using the meridian from zone 50 as a nonstandard meridian
         """
         tide=0
@@ -6821 +6828 @@
         urs2sts(base_name,
                 basename_out=base_name, 
-                zone=50,
-                mean_stage=tide,verbose=False)
+                central_meridian=123,
+                mean_stage=tide,
+                verbose=False)
 
         # now I want to check the sts file ...
@@ -6845 +6853 @@
         #Using the non standard projection (50) 
         for i in range(4):
-            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1], zone=50) 
+            zone, e, n = redfearn(lat_long_points[i][0],
+                                  lat_long_points[i][1],
+                                  central_meridian=123)
             assert num.allclose([x[i],y[i]], [e,n])
-            assert zone==geo_reference.zone
-
+            assert zone==-1
+
+            
     def test_urs2sts_nonstandard_projection_reverse(self):
         """
@@ -6907 +6918 @@
         #Using the non standard projection (50) 
         for i in range(4):
-            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1], zone=50) 
+            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1],
+                                  zone=50) 
             assert num.allclose([x[i],y[i]], [e,n])
             assert zone==geo_reference.zone
@@ -7191 +7203 @@
                 msg = 'urs stage is not equal to sts stage for for source %i and station %i' %(source_number,j)
                 assert num.allclose(urs_stage[index_start_urs_z:index_end_urs_z],
-                                sts_stage[index_start_stage:index_end_stage], 
-                                rtol=1.0e-6, atol=1.0e-5 ), msg                                
+                                    sts_stage[index_start_stage:index_end_stage], 
+                                    rtol=1.0e-6, atol=1.0e-5 ), msg                                
                                 
                 # check that urs e velocity and sts xmomentum are the same
@@ -7205 +7217 @@
                 msg = 'urs n velocity is not equivalent to sts y momentum for source %i and station %i' %(source_number,j)
                 assert num.allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]),
-                                sts_ymom[index_start_y:index_end_y], 
+                                -sts_ymom[index_start_y:index_end_y], 
                                 rtol=1.0e-5, atol=1.0e-4 ), msg
                                                 
@@ -10596 +10608 @@
         ymomentum = fid.variables['ymomentum'][:]        
 
+        
+        
         for i in range(stage.shape[0]):
             h = stage[i]-z # depth vector at time step i
@@ -10646 +10660 @@
 
         # 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#
+        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)

branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

@@ -2 +2 @@
 
 import unittest, os
+import os.path
 from math import pi, sqrt
 import tempfile
@@ -14 +15 @@
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
 
+from anuga.utilities.system_tools import get_pathname_from_package
 from shallow_water_domain import *
 
@@ -720 +722 @@
         # Check that flux seen from other triangles is inverse
         (ql, qr) = (qr, ql)
-        #tmp = qr
-        #qr = ql
-        #ql = tmp
         normal = domain.get_normal(2, 2)
         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
@@ -738 +737 @@
         # Check that flux seen from other triangles is inverse
         (ql, qr) = (qr, ql)
-        #tmp = qr
-        #qr = ql
-        #ql = tmp
         normal = domain.get_normal(3, 0)
         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
@@ -756 +752 @@
         # Check that flux seen from other triangles is inverse
         (ql, qr) = (qr, ql)
-        #tmp = qr
-        #qr=ql
-        #ql=tmp
         normal = domain.get_normal(0, 1)
         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
@@ -2472 +2465 @@
         R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
 
-        assert num.allclose(R.exchange_area, 1)
+        assert num.allclose(R.exchange_area, 2)
 
         domain.forcing_terms.append(R)
@@ -2508 +2501 @@
         # restricted to a polygon enclosing triangle #1 (bce)
         domain.forcing_terms = []
-        R = Rainfall(domain, rate=lambda t: 3*t + 7,
-                     polygon=[[1,1], [2,1], [2,2], [1,2]])
-
-        assert num.allclose(R.exchange_area, 1)
-
+        R = Rainfall(domain,
+                     rate=lambda t: 3*t + 7,
+                     polygon = [[1,1], [2,1], [2,2], [1,2]])
+
+        assert num.allclose(R.exchange_area, 2)
+        
         domain.forcing_terms.append(R)
 
@@ -2551 +2545 @@
         # restricted to a polygon enclosing triangle #1 (bce)
         domain.forcing_terms = []
-        R = Rainfall(domain, rate=lambda t: 3*t + 7,
-                     polygon=rainfall_poly)
-
-        assert num.allclose(R.exchange_area, 1)
-
+        R = Rainfall(domain,
+                     rate=lambda t: 3*t + 7,
+                     polygon=rainfall_poly)                     
+
+        assert num.allclose(R.exchange_area, 2)
+        
         domain.forcing_terms.append(R)
 
@@ -2616 +2611 @@
         # restricted to a polygon enclosing triangle #1 (bce)
         domain.forcing_terms = []
-        R = Rainfall(domain, rate=lambda t: 3*t + 7, polygon=rainfall_poly)
-
-        assert num.allclose(R.exchange_area, 1)
-
+        R = Rainfall(domain,
+                     rate=lambda t: 3*t + 7,
+                     polygon=rainfall_poly)
+
+        assert num.allclose(R.exchange_area, 2)
+        
         domain.forcing_terms.append(R)
 
@@ -2680 +2677 @@
 
         domain.forcing_terms = []
-        R = Rainfall(domain, rate=main_rate,
-                     polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0)
-
-        assert num.allclose(R.exchange_area, 1)
-
+        R = Rainfall(domain,
+                     rate=main_rate,
+                     polygon = [[1,1], [2,1], [2,2], [1,2]],
+                     default_rate=5.0)
+
+        assert num.allclose(R.exchange_area, 2)
+        
         domain.forcing_terms.append(R)
 
@@ -2746 +2745 @@
 
         domain.forcing_terms = []
-        R = Rainfall(domain, rate=main_rate,
-                     polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0)
-
-        assert num.allclose(R.exchange_area, 1)
-
+        R = Rainfall(domain,
+                     rate=main_rate,
+                     polygon=[[1,1], [2,1], [2,2], [1,2]],
+                     default_rate=5.0)
+
+        assert num.allclose(R.exchange_area, 2)
+        
         domain.forcing_terms.append(R)
 
@@ -2785 +2786 @@
         # on a circle affecting triangles #0 and #1 (bac and bce)
         domain.forcing_terms = []
-        domain.forcing_terms.append(Inflow(domain, rate=2.0,
-                                           center=(1,1), radius=1))
-
+        
+        I = Inflow(domain, rate=2.0, center=(1,1), radius=1)
+        domain.forcing_terms.append(I)        
         domain.compute_forcing_terms()
 
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            2.0/pi)
-        assert num.allclose(domain.quantities['stage'].explicit_update[0],
-                            2.0/pi)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
+        
+        A = I.exchange_area
+        assert num.allclose(A, 4) # Two triangles        
+        
+        assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A)
+        assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A)
+        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)        
+
 
     def test_inflow_using_circle_function(self):
@@ -2823 +2827 @@
         # on a circle affecting triangles #0 and #1 (bac and bce)
         domain.forcing_terms = []
-        domain.forcing_terms.append(Inflow(domain, rate=lambda t: 2.,
-                                           center=(1,1), radius=1))
-
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            2.0/pi)
-        assert num.allclose(domain.quantities['stage'].explicit_update[0],
-                            2.0/pi)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
+        I = Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1)
+        domain.forcing_terms.append(I)
+        
+        domain.compute_forcing_terms()        
+
+        A = I.exchange_area
+        assert num.allclose(A, 4) # Two triangles
+        
+        assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A)
+        assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A)
+        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)        
+        
+
+
 
     def test_inflow_catch_too_few_triangles(self):
@@ -6233 +6241 @@
         from anuga.geospatial_data.geospatial_data import Geospatial_data
 
-        #--------------------------------------------------------------------
+        
+        # Get path where this test is run
+        path = get_pathname_from_package('anuga.shallow_water')        
+        
+        
+        #----------------------------------------------------------------------
         # Create domain
         #--------------------------------------------------------------------
@@ -6276 +6289 @@
         interior_regions = []
         for polygon in offending_regions:
-            interior_regions.append( [polygon, 100] )
-
-        meshname = 'offending_mesh.msh'
+            interior_regions.append( [polygon, 100] ) 
+
+        meshname = os.path.join(path, 'offending_mesh.msh')
         create_mesh_from_regions(bounding_polygon,
                                  boundary_tags={'south': [0], 'east': [1],
@@ -6290 +6303 @@
         domain = Domain(meshname, use_cache=False, verbose=verbose)
 
-        #--------------------------------------------------------------------
+        #----------------------------------------------------------------------
         # Fit data point to mesh
-        #--------------------------------------------------------------------
-        points_file = 'offending_point.pts'
-
-        # This is the offending point
+        #----------------------------------------------------------------------
+
+        points_file = os.path.join(path, 'offending_point.pts')
+
+        # Offending point
         G = Geospatial_data(data_points=[[306953.344, 6194461.5]],
                             attributes=[1])
         G.export_points_file(points_file)
-
-        domain.set_quantity('elevation', filename=points_file, use_cache=False,
-                            verbose=verbose, alpha=0.01)
+        
+        try:
+            domain.set_quantity('elevation', filename=points_file,
+                                use_cache=False, verbose=verbose, alpha=0.01)
+        except RuntimeError, e:
+            msg = 'Test failed: %s' % str(e)
+            raise Exception, msg
+            # clean up in case raise fails
+            os.remove(meshname)
+            os.remove(points_file)
+        else:
+            os.remove(meshname)
+            os.remove(points_file)            
+        
+        #finally:
+            # Cleanup regardless
+            #FIXME(Ole): Finally does not work like this in python2.4 
+            #FIXME(Ole): Reinstate this when Python2.4 is out of the way
+            #FIXME(Ole): Python 2.6 apparently introduces something called 'with'            
+            #os.remove(meshname)
+            #os.remove(points_file)
+
+
+    def test_fitting_example_that_crashed_2(self):
+        """test_fitting_example_that_crashed_2
+        
+        This unit test has been derived from a real world example 
+        (the JJKelly study, by Petar Milevski).
+        
+        It shows a condition where set_quantity crashes due to AtA
+        not being built properly
+        
+        See ticket:314
+        """
+
+        verbose = False        
+
+        from anuga.shallow_water import Domain
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+        from anuga.geospatial_data import Geospatial_data
+        
+        # Get path where this test is run
+        path = get_pathname_from_package('anuga.shallow_water')        
+
+        meshname = os.path.join(path, 'test_mesh.msh')
+        W = 304180
+        S = 6185270
+        E = 307650
+        N = 6189040
+        maximum_triangle_area = 1000000
+
+        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
+
+        create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags={'south': [0], 
+                                                'east': [1], 
+                                                'north': [2], 
+                                                'west': [3]},
+                                 maximum_triangle_area=maximum_triangle_area,
+                                 filename=meshname,
+                                 use_cache=False,
+                                 verbose=verbose)
+
+        domain = Domain(meshname, use_cache=True, verbose=verbose)
+        
+        # Large test set revealed one problem
+        points_file = os.path.join(path, 'test_points_large.csv')
+        try:
+            domain.set_quantity('elevation', filename=points_file,
+                                use_cache=False, verbose=verbose)
+        except AssertionError, e:
+            msg = 'Test failed: %s' % str(e)
+            raise Exception, msg
+            # Cleanup in case this failed
+            os.remove(meshname)
+
+        # Small test set revealed another problem
+        points_file = os.path.join(path, 'test_points_small.csv')
+        try:    
+            domain.set_quantity('elevation', filename=points_file,
+                                use_cache=False, verbose=verbose)                            
+        except AssertionError, e:
+            msg = 'Test failed: %s' % str(e)
+            raise Exception, msg
+            # Cleanup in case this failed
+            os.remove(meshname)
+        else:
+            os.remove(meshname)
+
+            
+    def test_total_volume(self):        
+        """test_total_volume
+        
+        Test that total volume can be computed correctly
+        """            
+
+        #----------------------------------------------------------------------
+        # Import necessary modules
+        #----------------------------------------------------------------------
+        from anuga.abstract_2d_finite_volumes.mesh_factory \
+                import rectangular_cross
+        from anuga.shallow_water import Domain
+
+        #----------------------------------------------------------------------
+        # Setup computational domain
+        #----------------------------------------------------------------------
+
+        length = 100.
+        width  = 20.
+        dx = dy = 5       # Resolution: of grid on both axes
+        
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length,
+                                                       len2=width)
+        domain = Domain(points, vertices, boundary)   
+
+        #----------------------------------------------------------------------
+        # Simple flat bottom bathtub
+        #----------------------------------------------------------------------
+
+        d = 1.0
+        domain.set_quantity('elevation', 0.0)
+        domain.set_quantity('stage', d)
+        
+        assert num.allclose(domain.compute_total_volume(), length*width*d)
+
+        #----------------------------------------------------------------------
+        # Slope
+        #----------------------------------------------------------------------
+                
+        slope = 1.0/10          # RHS drops to -10m
+        def topography(x, y):
+            return -x * slope
+
+        domain.set_quantity('elevation', topography)
+        domain.set_quantity('stage', 0.0)       # Domain full
+        
+        V = domain.compute_total_volume()
+        assert num.allclose(V, length*width*10/2)
+
+        domain.set_quantity('stage', -5.0)      # Domain 'half' full
+        
+        # IMPORTANT: Adjust stage to match elevation
+        domain.distribute_to_vertices_and_edges()
+        
+        V = domain.compute_total_volume()
+        assert num.allclose(V, width*(length/2)*5.0/2)
+
+
+    def test_volumetric_balance_computation(self):
+        """test_volumetric_balance_computation
+        
+        Test that total in and out flows are computed correctly
+        FIXME(Ole): This test is more about looking at the printed report
+        """
+
+        verbose = False
+
+        #----------------------------------------------------------------------
+        # Import necessary modules
+        #----------------------------------------------------------------------
+
+        from anuga.abstract_2d_finite_volumes.mesh_factory \
+                import rectangular_cross
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+        from anuga.shallow_water.shallow_water_domain import Inflow
+        from anuga.shallow_water.data_manager \
+                import get_flow_through_cross_section
+
+        #----------------------------------------------------------------------
+        # Setup computational domain
+        #----------------------------------------------------------------------
+
+        finaltime = 300.0
+        length = 300.
+        width  = 20.
+        dx = dy = 5       # Resolution: of grid on both axes
+        
+        # Input parameters
+        uh = 1.0
+        vh = 0.0
+        d = 1.0
+        
+        # 20 m^3/s in the x direction across entire domain
+        ref_flow = uh*d*width
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length,
+                                                       len2=width)
+
+        domain = Domain(points, vertices, boundary)   
+        domain.set_name('Inflow_flowline_test')              # Output name
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+
+        slope = 0.0
+        def topography(x, y):
+            z=-x * slope
+            return z
+
+        # Use function for elevation
+        domain.set_quantity('elevation', topography)
+        domain.set_quantity('friction', 0.0)        # Constant friction
+                
+        domain.set_quantity('stage', expression='elevation')
+
+        #----------------------------------------------------------------------
+        # Setup boundary conditions
+        #----------------------------------------------------------------------
+
+        Br = Reflective_boundary(domain)      # Solid reflective wall
+                
+        # Constant flow into domain
+        # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s 
+        
+        Bi = Dirichlet_boundary([d, uh, vh]) 
+        Bo = Dirichlet_boundary([0, 0, 0])
+
+        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
+
+        #----------------------------------------------------------------------
+        # Evolve system through time
+        #----------------------------------------------------------------------
+
+        for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
+            S = domain.volumetric_balance_statistics()
+            if verbose :
+                print domain.timestepping_statistics()
+                print S
+
+
+    def Xtest_inflow_using_flowline(self):
+        """test_inflow_using_flowline
+
+        Test the ability of a flowline to match inflow above the flowline by
+        creating constant inflow onto a circle at the head of a 20m wide by 300m
+        long plane dipping at various slopes with a perpendicular flowline and
+        gauge downstream of the inflow and a 45 degree flowlines at 200m
+        downstream.
+        """
+
+        # FIXME(Ole): Move this and the following test to validate_all.py as
+        # they are rather time consuming
+        
+        verbose = True
+
+        #----------------------------------------------------------------------
+        # Import necessary modules
+        #----------------------------------------------------------------------
+
+        from anuga.abstract_2d_finite_volumes.mesh_factory \
+                import rectangular_cross
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+        from anuga.shallow_water.shallow_water_domain import Inflow
+        from anuga.shallow_water.data_manager \
+                import get_flow_through_cross_section
+        from anuga.abstract_2d_finite_volumes.util \
+                import sww2csv_gauges, csv2timeseries_graphs
+
+        #----------------------------------------------------------------------
+        # Setup computational domain
+        #----------------------------------------------------------------------
+
+        number_of_inflows = 2       # Number of inflows on top of each other
+        finaltime = 1000.0          # 6000.0
+
+        length = 300.
+        width  = 20.
+        dx = dy = 5                 # Resolution: of grid on both axes
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length,
+                                                       len2=width)
+
+        for mannings_n in [0.150, 0.07, 0.035]:    #[0.012, 0.035, 0.070, 0.150]:
+            for slope in [1.0/300, 1.0/150, 1.0/75]:
+                # Loop over a range of bedslopes representing sub to
+                # super critical flows 
+                if verbose:
+                    print
+                    print 'Slope:', slope, 'Mannings n:', mannings_n
+                domain = Domain(points, vertices, boundary)   
+                domain.set_name('Inflow_flowline_test')     # Output name
+
+                #--------------------------------------------------------------
+                # Setup initial conditions
+                #--------------------------------------------------------------
+
+                def topography(x, y):
+                    z = -x * slope
+                    return z
+
+                # Use function for elevation
+                domain.set_quantity('elevation', topography)
+                # Constant friction of conc surface   
+                domain.set_quantity('friction', mannings_n)
+                # Dry initial condition
+                domain.set_quantity('stage', expression='elevation')
+
+                #--------------------------------------------------------------
+                # Setup boundary conditions
+                #--------------------------------------------------------------
+
+                Br = Reflective_boundary(domain)      # Solid reflective wall
+                # Outflow stage at -0.9m d=0.1m
+                Bo = Dirichlet_boundary([-100, 0, 0])
+
+                domain.set_boundary({'left': Br, 'right': Bo,
+                                     'top': Br,  'bottom': Br})
+
+                #--------------------------------------------------------------
+                # Setup Inflow
+                #--------------------------------------------------------------
+
+                # Fixed Flowrate onto Area 
+                fixed_inflow = Inflow(domain, center=(10.0, 10.0),
+                                      radius=5.00, rate=10.00)   
+
+                # Stack this flow
+                for i in range(number_of_inflows):
+                    domain.forcing_terms.append(fixed_inflow)
+                
+                ref_flow = fixed_inflow.rate*number_of_inflows
+
+                #--------------------------------------------------------------
+                # Evolve system through time
+                #--------------------------------------------------------------
+
+                for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
+                    pass
+                    #if verbose :
+                    #    print domain.timestepping_statistics()
+
+                if verbose:
+                    print domain.volumetric_balance_statistics()                                                            
+                #--------------------------------------------------------------
+                # Compute flow thru flowlines ds of inflow
+                #--------------------------------------------------------------
+                    
+                # Square on flowline at 200m
+                q = domain.get_flow_through_cross_section([[200.0,  0.0],
+                                                           [200.0, 20.0]])
+                if verbose:
+                    print ('90 degree flowline: ANUGA = %f, Ref = %f'
+                           % (q, ref_flow))
+
+                msg = ('Predicted flow was %f, should have been %f'
+                       % (q, ref_flow))
+                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
+
+                           
+                # 45 degree flowline at 200m
+                q = domain.get_flow_through_cross_section([[200.0,  0.0],
+                                                           [220.0, 20.0]])
+                if verbose:
+                    print ('45 degree flowline: ANUGA = %f, Ref = %f'
+                           % (q, ref_flow))
+                    
+                msg = ('Predicted flow was %f, should have been %f'
+                       % (q, ref_flow))
+                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
+
+                # Stage recorder (gauge) in middle of plane at 200m
+                x = 200.0
+                y = 10.00
+                w = domain.get_quantity('stage').\
+                        get_values(interpolation_points=[[x, y]])[0]
+                z = domain.get_quantity('elevation').\
+                        get_values(interpolation_points=[[x, y]])[0]
+                domain_depth = w-z
+
+                # Compute normal depth at gauge location using Manning equation
+                # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
+                normal_depth = (ref_flow*mannings_n/(slope**0.5*width))**0.6
+                msg = ('Predicted depth of flow was %f, should have been %f'
+                       % (domain_depth, normal_depth))
+                if verbose:
+                    diff = abs(domain_depth-normal_depth)
+                    print ('Depth: ANUGA = %f, Mannings = %f, E=%f'
+                           % (domain_depth, normal_depth, diff/domain_depth))
+                    assert diff < 0.1
+                    
+                if slope >= 1.0/100:
+                    # Really super critical flow is not as stable.
+                    #assert num.allclose(domain_depth, normal_depth,
+                    #                    rtol=1.0e-1), msg
+                    pass
+                else:
+                    pass
+                    #assert num.allclose(domain_depth, normal_depth,
+                    #                    rtol=1.0e-2), msg
+
+
+    def Xtest_friction_dependent_flow_using_flowline(self):
+        """test_friction_dependent_flow_using_flowline
+        
+        Test the internal flow (using flowline) as a function of
+        different values of Mannings n and different slopes.
+        
+        Flow is applied in the form of boundary conditions with fixed momentum.
+        """
+
+        verbose = True
+
+        #----------------------------------------------------------------------
+        # Import necessary modules
+        #----------------------------------------------------------------------
+
+        from anuga.abstract_2d_finite_volumes.mesh_factory \
+                import rectangular_cross
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+        from anuga.shallow_water.shallow_water_domain import Inflow
+        from anuga.shallow_water.data_manager \
+                import get_flow_through_cross_section
+        from anuga.abstract_2d_finite_volumes.util \
+                import sww2csv_gauges, csv2timeseries_graphs
+
+
+        #----------------------------------------------------------------------
+        # Setup computational domain
+        #----------------------------------------------------------------------
+
+        finaltime = 1000.0
+
+        length = 300.
+        width  = 20.
+        dx = dy = 5       # Resolution: of grid on both axes
+        
+        # Input parameters
+        uh = 1.0
+        vh = 0.0
+        d = 1.0
+        
+        ref_flow = uh*d*width # 20 m^3/s in the x direction across entire domain
+
+        points, vertices, boundary = rectangular_cross(int(length/dx),
+                                                       int(width/dy),
+                                                       len1=length,
+                                                       len2=width)
+
+        for mannings_n in [0.035]:          #[0.0, 0.012, 0.035]:
+            for slope in [1.0/300]:         #[0.0, 1.0/300, 1.0/150]:
+                # Loop over a range of bedslopes representing
+                # sub to super critical flows 
+                if verbose:
+                    print
+                    print 'Slope:', slope, 'Mannings n:', mannings_n
+                domain = Domain(points, vertices, boundary)   
+                domain.set_name('Inflow_flowline_test')     # Output name
+
+                #--------------------------------------------------------------
+                # Setup initial conditions
+                #--------------------------------------------------------------
+
+                def topography(x, y):
+                    z = -x * slope
+                    return z
+
+                # Use function for elevation
+                domain.set_quantity('elevation', topography)
+                # Constant friction
+                domain.set_quantity('friction', mannings_n)
+                
+                #domain.set_quantity('stage', expression='elevation')
+                     
+                # Set initial flow as depth=1m, uh=1.0 m/s, vh = 0.0
+                # making it 20 m^3/s across entire domain 
+                domain.set_quantity('stage', expression='elevation + %f' % d)
+                domain.set_quantity('xmomentum', uh)
+                domain.set_quantity('ymomentum', vh)                
+
+                #--------------------------------------------------------------
+                # Setup boundary conditions
+                #--------------------------------------------------------------
+
+                Br = Reflective_boundary(domain)      # Solid reflective wall
+                
+                # Constant flow in and out of domain
+                # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s 
+                # across boundaries
+                Bi = Dirichlet_boundary([d, uh, vh]) 
+                Bo = Dirichlet_boundary([-length*slope+d, uh, vh])
+                #Bo = Dirichlet_boundary([-100, 0, 0])
+
+                domain.set_boundary({'left': Bi, 'right': Bo,
+                                     'top': Br,  'bottom': Br})
+
+                #--------------------------------------------------------------
+                # Evolve system through time
+                #--------------------------------------------------------------
+
+                for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
+                    if verbose :
+                        print domain.timestepping_statistics()
+                        print domain.volumetric_balance_statistics()
+
+                # 90 degree flowline at 200m
+                q = domain.get_flow_through_cross_section([[200.0,  0.0],
+                                                           [200.0, 20.0]])
+                msg = ('Predicted flow was %f, should have been %f'
+                       % (q, ref_flow))
+                if verbose:
+                    print ('90 degree flowline: ANUGA = %f, Ref = %f'
+                           % (q, ref_flow))
+
+                # 45 degree flowline at 200m
+                q = domain.get_flow_through_cross_section([[200.0,  0.0],
+                                                           [220.0, 20.0]])
+                msg = ('Predicted flow was %f, should have been %f'
+                       % (q, ref_flow))
+                if verbose:
+                    print ('45 degree flowline: ANUGA = %f, Ref = %f'
+                           % (q, ref_flow))
+
+                # Stage recorder (gauge) in middle of plane at 200m
+                x = 200.0
+                y = 10.00                
+                w = domain.get_quantity('stage').\
+                        get_values(interpolation_points=[[x, y]])[0]
+                z = domain.get_quantity('elevation').\
+                        get_values(interpolation_points=[[x, y]])[0]
+                domain_depth = w-z
+
+                xmom = domain.get_quantity('xmomentum').\
+                        get_values(interpolation_points=[[x, y]])[0]
+                ymom = domain.get_quantity('ymomentum').\
+                        get_values(interpolation_points=[[x, y]])[0]            
+                if verbose:                    
+                    print ('At interpolation point (h, uh, vh): ',
+                           domain_depth, xmom, ymom)
+                    print 'uh * d * width = ', xmom*domain_depth*width
+                                
+                if slope > 0.0:
+                    # Compute normal depth at gauge location using Manning eqn
+                    # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
+                    normal_depth = (ref_flow*mannings_n/(slope**0.5*width))**0.6
+                    if verbose:
+                        print ('Depth: ANUGA = %f, Mannings = %f'
+                               % (domain_depth, normal_depth))
+
+#################################################################################
 
 if __name__ == "__main__":

branches/numpy/anuga/utilities/log.py

@@ -7 +7 @@
 Use it this way:
     import anuga.utilities.log as log
+    log.console_logging_level = log.DEBUG
     log.debug('A message at DEBUG level')
     log.info('Another message, INFO level')
@@ -17 +18 @@
 Until the first call to log() the user is free to play with the module data
 to configure the logging.
+Note that this module uses features of the logging package that were introduced
+in python2.5.  If running on earlier versions, these features are disables:
+    . Module name + line number
 '''
 
 import sys
-import os.path
+import os
+import sys
 import traceback
 import logging
@@ -27 +32 @@
 ################################################################################
 # Module variables - only one copy of these, ever.
+#
+# The console logging level is set to a high level, like CRITICAL.  The logfile
+# logging is set lower, between DEBUG and CRITICAL.  The idea is to log least to
+# the console, but ensure that everything that goes to the console *will* also
+# appear in the log file.  There is code to ensure log <= console levels.
 ################################################################################
 
@@ -33 +43 @@
 
 # logging level for the console
-console_logging_level = logging.INFO
+console_logging_level = logging.CRITICAL
 
 # logging level for the logfile
-log_logging_level = logging.DEBUG
-
-# The name of the file to log to.
-log_filename = './anuga.log'
+log_logging_level = logging.INFO
+
+# The default name of the file to log to.
+log_filename = os.path.join('.', 'anuga.log')
+
+# set module variables so users don't have to do 'import logging'.
+CRITICAL = logging.CRITICAL
+ERROR = logging.ERROR
+WARNING = logging.WARNING
+INFO = logging.INFO
+DEBUG = logging.DEBUG
+NOTSET = logging.NOTSET
 
 
@@ -58 +76 @@
     '''
 
-    global _setup
+    global _setup, log_logging_level
+
+    # get running python version for later
+    (version_major, version_minor, _, _, _) = sys.version_info
 
     # have we been setup?
     if not _setup:
+        # sanity check the logging levels, require console >= file
+        if console_logging_level < log_logging_level:
+            log_logging_level = console_logging_level
+
         # setup the file logging system
-        fmt = '%(asctime)s %(levelname)-8s %(mname)25s:%(lnum)-4d|%(message)s'
+        if version_major >= 2 and version_minor >= 5:
+            fmt = '%(asctime)s %(levelname)-8s %(mname)25s:%(lnum)-4d|%(message)s'
+        else:
+            fmt = '%(asctime)s %(levelname)-8s|%(message)s'
         logging.basicConfig(level=log_logging_level, format=fmt,
                             filename=log_filename, filemode='w')
@@ -80 +108 @@
                         logging.getLevelName(log_logging_level),
                         logging.getLevelName(console_logging_level)))
-        logging.log(logging.CRITICAL, start_msg,
-                    extra={'mname': __name__, 'lnum': 0})
+        if version_major >= 2 and version_minor >= 5:
+            logging.log(logging.CRITICAL, start_msg,
+                        extra={'mname': __name__, 'lnum': 0})
+        else:
+            logging.log(logging.CRITICAL, start_msg)
 
         # mark module as *setup*
@@ -89 +120 @@
     frames = traceback.extract_stack()
     frames.reverse()
-    for (mname, lnum, _, _) in frames:
-        mname = os.path.basename(mname).rsplit('.', 1)[0]
+    for (fpath, lnum, mname, _) in frames:
+        fname = os.path.basename(mname).rsplit('.', 1)[0]
         if mname != __name__:
             break
 
-    logging.log(level, msg, extra={'mname': mname, 'lnum': lnum})
+    if version_major >= 2 and version_minor >= 5:
+        logging.log(level, msg, extra={'mname': fname, 'lnum': lnum})
+    else:
+        logging.log(level, msg)
 
 ################################################################################
@@ -130 +164 @@
     log(logging.CRITICAL, msg)
 
+def resource_usage(level=logging.CRITICAL):
+    '''Log resource usage at given log level.'''
+
+    if sys.platform != 'win32':
+        _proc_status = '/proc/%d/status' % os.getpid()
+        _scale = {'KB': 1024.0, 'MB': 1024.0*1024.0, 'GB': 1024.0*1024.0*1024.0,
+                  'kB': 1024.0, 'mB': 1024.0*1024.0, 'gB': 1024.0*1024.0*1024.0}
+
+        def _VmB(VmKey):
+            '''Get number of virtual bytes used.'''
+
+            # get pseudo file /proc/<pid>/status
+            try:
+                t = open(_proc_status)
+                v = t.read()
+                t.close()
+            except IOError:
+                return 0.0
+
+            # get VmKey line, eg: 'VmRSS: 999 kB\n ...
+            i = v.index(VmKey)
+            v = v[i:].split(None, 3)
+            if len(v) < 3:
+                return 0.0
+
+            # convert Vm value to bytes
+            return float(v[1]) * _scale[v[2]]
+
+        def memory(since=0.0):
+            '''Get virtual memory usage in bytes.'''
+
+            return _VmB('VmSize:') - since
+
+        def resident(since=0.0):
+            '''Get resident memory usage in bytes.'''
+
+            return _VmB('VmRSS:') - since
+
+        def stacksize(since=0.0):
+            '''Get stack size in bytes.'''
+
+            return _VmB('VmStk:') - since
+
+        msg = ('Resource usage: memory=%.1f resident=%.1f stacksize=%.1f'
+               % (memory()/_scale['GB'], resident()/_scale['GB'],
+                  stacksize()/_scale['GB']))
+        log(level, msg)
+    else:
+        pass

branches/numpy/anuga/utilities/polygon.py

@@ -1070 +1070 @@
     class Found(exceptions.Exception): pass
 
+    polygon = ensure_numeric(polygon)
+    
     point_in = False
     while not point_in:
@@ -1119 +1121 @@
     # TO DO check if any of the regions fall inside one another
 
-    print '--------------------------------------------------------------------'
-    print ('Polygon   Max triangle area (m^2)   Total area (km^2)   '
-           'Estimated #triangles')
-    print '--------------------------------------------------------------------'
-
+    print '----------------------------------------------------------------------------'
+    print 'Polygon   Max triangle area (m^2)   Total area (km^2)   Estimated #triangles'
+    print '----------------------------------------------------------------------------'    
+        
     no_triangles = 0.0
     area = polygon_area(bounding_poly)

branches/numpy/anuga/utilities/system_tools.py

@@ -7 +7 @@
 import os
 import string
+import urllib
+import urllib2
+import getpass
+import tarfile
+import md5
+
 
 def log_to_file(filename, s, verbose=False):
@@ -49 +55 @@
     """
 
-    # Create dummy info 
-    #info = 'Revision: Version info could not be obtained.'
-    #info += 'A command line version of svn must be availbable '
-    #info += 'on the system PATH, access to the subversion '
-    #info += 'repository is necessary and the output must '
-    #info += 'contain a line starting with "Revision:"'
-    
-
-    #FIXME (Ole): Change this so that svn info is attempted first.
-    # If that fails, try to read a stored file with that same info (this would be created by e.g. the release script). Failing that, throw an exception.
-
-    #FIXME (Ole): Move this and store_version_info to utilities
-
-
+    # FIXME (Ole): Change this so that svn info is attempted first.
+    # If that fails, try to read a stored file with that same info (this would
+    # be created by e.g. the release script). Failing that, throw an exception.
+    # FIXME (Ole): Move this and store_version_info to utilities
     try:
         from anuga.stored_version_info import version_info
     except:
-        msg = 'No version info stored and command "svn" is not '
-        msg += 'recognised on the system PATH.\n\n'
-        msg += 'If ANUGA has been installed from a distribution e.g. as '
-        msg += 'obtained from SourceForge,\n'
-        msg += 'the version info should be '
-        msg += 'available in the automatically generated file '
-        msg += 'stored_version_info.py\n'
-        msg += 'in the anuga root directory.\n'
-        msg += 'If run from a Subversion sandpit, '
-        msg += 'ANUGA will try to obtain the version info '
-        msg += 'by using the command: "svn info".\n'
-        msg += 'In this case, make sure svn is accessible on the system path. '
-        msg += 'Simply aliasing svn to the binary will not work. '
-        msg += 'Good luck!'
+        msg = '''
+No version info stored and command 'svn' is not recognised on the system PATH.
+
+If ANUGA has been installed from a distribution e.g. as obtained from
+SourceForge, the version info should be available in the automatically
+generated file 'stored_version_info.py' in the anuga root directory.
+
+If run from a Subversion sandpit, ANUGA will try to obtain the version info
+by using the command 'svn info'.  In this case, make sure 'svn' is accessible
+on the system PATH.  Simply aliasing 'svn' to the binary will not work.
+
+If you are using Windows, you have to install the file svn.exe which can be
+obtained from http://www.collab.net/downloads/subversion.
+
+Good luck!
+'''
 
         # No file available - try using Subversion
@@ -91 +89 @@
             else:
                 fid = os.popen('svn info 2>/dev/null')
-        
         except:
             raise Exception(msg)
         else:
-            #print 'Got version from svn'            
             version_info = fid.read()
-            
-            if version_info == '':
-                raise Exception(msg)    
+
+        if version_info == '':
+            raise Exception, msg
     else:
         pass
         #print 'Got version from file'
-
             
     for line in version_info.split('\n'):
@@ -110 +105 @@
 
     fields = line.split(':')
-    msg = 'Keyword "Revision" was not found anywhere in text: %s' %version_info
+    msg = 'Keyword "Revision" was not found anywhere in text: %s' % version_info
     assert fields[0].startswith('Revision'), msg            
 
@@ -118 +113 @@
         msg = 'Revision number must be an integer. I got %s' %fields[1]
         msg += 'Check that the command svn is on the system path' 
-        raise Exception(msg)                
-        
+        raise Exception, msg
+
     return revision_number
 
@@ -181 +176 @@
         if verbose is True:
             print 'Version info stored to %s' %filename
+
 
 def safe_crc(string):
@@ -319 +315 @@
 
 
+##
+# @brief Get list of variable names in a python expression string.
+# @param source A string containing a python expression.
+# @return A list of variable name strings.
+# @note Throws SyntaxError exception if not a valid expression.
+def get_vars_in_expression(source):
+    '''Get list of variable names in a python expression.'''
+
+    import compiler
+    from compiler.ast import Node
+
+    ##
+    # @brief Internal recursive function.
+    # @param node An AST parse Node.
+    # @param var_list Input list of variables.
+    # @return An updated list of variables.
+    def get_vars_body(node, var_list=[]):
+        if isinstance(node, Node):
+            if node.__class__.__name__ == 'Name':
+                for child in node.getChildren():
+                    if child not in var_list:
+                        var_list.append(child)
+            for child in node.getChildren():
+                if isinstance(child, Node):
+                    for child in node.getChildren():
+                        var_list = get_vars_body(child, var_list)
+                    break
+
+        return var_list
+
+    return get_vars_body(compiler.parse(source))
+
+
+##
+# @brief Get a file from the web.
+# @param file_url URL of the file to fetch.
+# @param file_name Path to file to create in the filesystem.
+# @param auth Auth tuple (httpproxy, proxyuser, proxypass).
+# @param blocksize Read file in this block size.
+# @return 'auth' tuple for subsequent calls, if successful.
+# @note If 'auth' not supplied, will prompt user.
+# @note Will try using environment variable HTTP_PROXY for proxy server.
+# @note Will try using environment variable PROXY_USERNAME for proxy username.
+# @note Will try using environment variable PROXY_PASSWORD for proxy password.
+def get_web_file(file_url, file_name, auth=None, blocksize=1024*1024):
+    '''Get a file from the web.
+
+    file_url:  The URL of the file to get
+    file_name: Local path to save loaded file in
+    auth:      A tuple (httpproxy, proxyuser, proxypass)
+    blocksize: Block size of file reads
+    
+    Will try simple load through urllib first.  Drop down to urllib2
+    if there is a proxy and it requires authentication.
+
+    Environment variable HTTP_PROXY can be used to supply proxy information.
+    PROXY_USERNAME is used to supply the authentication username.
+    PROXY_PASSWORD supplies the password, if you dare!
+    '''
+
+    # Simple fetch, if fails, check for proxy error
+    try:
+        urllib.urlretrieve(file_url, file_name)
+        return None     # no proxy, no auth required
+    except IOError, e:
+        if e[1] != 407:
+            raise       # raise error if *not* proxy auth error
+
+    # We get here if there was a proxy error, get file through the proxy
+    # unpack auth info
+    try:
+        (httpproxy, proxyuser, proxypass) = auth
+    except:
+        (httpproxy, proxyuser, proxypass) = (None, None, None)
+
+    # fill in any gaps from the environment
+    if httpproxy is None:
+        httpproxy = os.getenv('HTTP_PROXY')
+    if proxyuser is None:
+        proxyuser = os.getenv('PROXY_USERNAME')
+    if proxypass is None:
+        proxypass = os.getenv('PROXY_PASSWORD')
+
+    # Get auth info from user if still not supplied
+    if httpproxy is None or proxyuser is None or proxypass is None:
+        print '-'*80
+        print ('You need to supply proxy authentication information.  '
+               'Use environment variables\n'
+               'HTTP_PROXY, PROXY_USERNAME and PROXY_PASSWORD to bypass '
+               'entry here:')
+        if httpproxy is None:
+            httpproxy = raw_input('  proxy server: ')
+        if proxyuser is None:
+            proxyuser = raw_input('proxy username: ') 
+        if proxypass is None:
+            proxypass = getpass.getpass('proxy password: ')
+        print '-'*80
+
+    # the proxy URL cannot start with 'http://', we add that later
+    httpproxy = httpproxy.lower()
+    if httpproxy.startswith('http://'):
+        httpproxy = httpproxy.replace('http://', '', 1)
+
+    # open remote file
+    proxy = urllib2.ProxyHandler({'http': 'http://' + proxyuser
+                                              + ':' + proxypass
+                                              + '@' + httpproxy})
+    authinfo = urllib2.HTTPBasicAuthHandler()
+    opener = urllib2.build_opener(proxy, authinfo, urllib2.HTTPHandler)
+    urllib2.install_opener(opener)
+    webget = urllib2.urlopen(file_url)
+
+    # transfer file to local filesystem
+    fd = open(file_name, 'wb')
+    while True:
+        data = webget.read(blocksize)
+        if len(data) == 0:
+            break
+        fd.write(data)
+    fd.close
+    webget.close()
+
+    # return successful auth info
+    return (httpproxy, proxyuser, proxypass)
+
+
+##
+# @brief Tar a file (or directory) into a tarfile.
+# @param files A list of files (or directories) to tar.
+# @param tarfile The created tarfile name.
+# @note We use gzip compression.
+def tar_file(files, tarname):
+    '''Compress a file or directory into a tar file.'''
+
+    o = tarfile.open(tarname, 'w:gz')
+    for file in files:
+        o.add(file)
+    o.close()
+
+
+##
+# @brief Untar a file into an optional target directory.
+# @param tarname Name of the file to untar.
+# @param target_dir Directory to untar into.
+def untar_file(tarname, target_dir='.'):
+    '''Uncompress a tar file.'''
+
+    o = tarfile.open(tarname, 'r:gz')
+    members = o.getmembers()
+    for member in members:
+        o.extract(member, target_dir)
+    o.close()
+
+
+##
+# @brief Return a hex digest string for a given file.
+# @param filename Path to the file of interest.
+# @param blocksize Size of data blocks to read.
+# @return A hex digest string (16 bytes).
+# @note Uses MD5 digest.
+def get_file_hexdigest(filename, blocksize=1024*1024*10):
+    '''Get a hex digest of a file.'''
+    
+    m = md5.new()
+    fd = open(filename, 'r')
+            
+    while True:
+        data = fd.read(blocksize)
+        if len(data) == 0:
+            break
+        m.update(data)
+                                                                
+    fd.close()
+    return m.hexdigest()
+
+    fd = open(filename, 'r')
+
+
+##
+# @brief Create a file containing a hexdigest string of a data file.
+# @param data_file Path to the file to get the hexdigest from.
+# @param digest_file Path to hexdigest file to create.
+# @note Uses MD5 digest.
+def make_digest_file(data_file, digest_file):
+    '''Create a file containing the hex digest string of a data file.'''
+    
+    hexdigest = get_file_hexdigest(data_file)
+    fd = open(digest_file, 'w')
+    fd.write(hexdigest)
+    fd.close()
+
+
+##
+# @brief Function to return the length of a file.
+# @param in_file Path to file to get length of.
+# @return Number of lines in file.
+# @note Doesn't count '\n' characters.
+# @note Zero byte file, returns 0.
+# @note No \n in file at all, but >0 chars, returns 1.
+def file_length(in_file):
+    '''Function to return the length of a file.'''
+
+    fid = open(in_file)
+    data = fid.readlines()
+    fid.close()
+    return len(data)
+
+

branches/numpy/anuga/utilities/test_polygon.py

@@ -1381 +1381 @@
             assert status == 1
 
-            msg = ('Orientation of line shouldn not matter.\n'
+            msg = ('Orientation of line should not matter.\n'
                    'However, segment [%f,%f], [%f, %f]' %
                    (x, 0, common_end_point[0], common_end_point[1]))

branches/numpy/anuga/utilities/test_system_tools.py

@@ -4 +4 @@
 import unittest
 import numpy as num
+import random
+import tempfile
 import zlib
 import os
@@ -85 +87 @@
             pass
         else:
-            test_array = num.array([[7.0, 3.14], [-31.333, 0.0]],
-                                   dtype=num.float)
+            test_array = num.array([[7.0, 3.14], [-31.333, 0.0]])
 
             # First file
@@ -305 +306 @@
         os.remove(FILENAME)
 
-#        print ('test_string_to_char: str_list=%s, new_str_list=%s'
-                #               % (str(str_list), str(new_str_list)))
+
+    def test_get_vars_in_expression(self):
+        '''Test the 'get vars from expression' code.'''
+
+        def test_it(source, expected):
+            result = get_vars_in_expression(source)
+            result.sort()
+            expected.sort()
+            msg = ("Source: '%s'\nResult: %s\nExpected: %s"
+                   % (source, str(result), str(expected)))
+            self.failUnlessEqual(result, expected, msg)
+                
+        source = 'fred'
+        expected = ['fred']
+        test_it(source, expected)
+
+        source = 'tom + dick'
+        expected = ['tom', 'dick']
+        test_it(source, expected)
+
+        source = 'tom * (dick + harry)'
+        expected = ['tom', 'dick', 'harry']
+        test_it(source, expected)
+
+        source = 'tom + dick**0.5 / (harry - tom)'
+        expected = ['tom', 'dick', 'harry']
+        test_it(source, expected)
+
+
+    def test_tar_untar_files(self):
+        '''Test that tarring & untarring files is OK.'''
+
+        num_lines = 100
+        line_size = 100
+
+        # these test files must exist in the current directory
+        # create them with random data
+        files = ('alpha', 'beta', 'gamma')
+        for file in files:
+            fd = open(file, 'w')
+            line = ''
+            for i in range(num_lines):
+                for j in range(line_size):
+                    line += chr(random.randint(ord('A'), ord('Z')))
+                line += '\n'
+                fd.write(line)
+            fd.close()
+
+        # name of tar file and test (temp) directory
+        tar_filename = 'test.tgz'
+        tmp_dir = tempfile.mkdtemp()
+
+        # tar and untar the test files into a temporary directory
+        tar_file(files, tar_filename)
+        untar_file(tar_filename, tmp_dir)
+
+        # see if original files and untarred ones are the same
+        for file in files:
+            fd = open(file, 'r')
+            orig = fd.readlines()
+            fd.close()
+
+            fd = open(os.path.join(tmp_dir, file), 'r')
+            copy = fd.readlines()
+            fd.close()
+
+            msg = "Original file %s isn't the same as untarred copy?" % file
+            self.failUnless(orig == copy, msg)
+
+        # clean up
+        for file in files:
+            os.remove(file)
+        os.remove(tar_filename)
+
+
+    def test_file_digest(self):
+        '''Test that file digest functions give 'correct' answer.
+        
+        Not a good test as we get 'expected_digest' from a digest file,
+        but *does* alert us if the digest algorithm ever changes.
+        '''
+
+        # we expect this digest string from the data file
+        expected_digest = '831a1dde6edd365ec4163a47871fa21b'
+
+        # prepare test directory and filenames
+        tmp_dir = tempfile.mkdtemp()
+        data_file = os.path.join(tmp_dir, 'test.data')
+        digest_file = os.path.join(tmp_dir, 'test.digest')
+
+        # create the data file
+        data_line = 'The quick brown fox jumps over the lazy dog. 0123456789\n'
+        fd = open(data_file, 'w')
+        for line in range(100):
+            fd.write(data_line)
+        fd.close()
+
+        # create the digest file
+        make_digest_file(data_file, digest_file)
+
+        # get digest string for the data file
+        digest = get_file_hexdigest(data_file)
+
+        # check that digest is as expected, string
+        msg = ("Digest string wrong, got '%s', expected '%s'"
+               % (digest, expected_digest))
+        self.failUnless(expected_digest == digest, msg)
+
+        # check that digest is as expected, file
+        msg = ("Digest file wrong, got '%s', expected '%s'"
+               % (digest, expected_digest))
+        fd = open(digest_file, 'r')
+        digest = fd.readline()
+        fd.close()
+        self.failUnless(expected_digest == digest, msg)
+
+
+    def test_file_length_function(self):
+        '''Test that file_length() give 'correct' answer.'''
+
+        # prepare test directory and filenames
+        tmp_dir = tempfile.mkdtemp()
+        test_file1 = os.path.join(tmp_dir, 'test.file1')
+        test_file2 = os.path.join(tmp_dir, 'test.file2')
+        test_file3 = os.path.join(tmp_dir, 'test.file3')
+        test_file4 = os.path.join(tmp_dir, 'test.file4')
+
+        # create files of known length
+        fd = open(test_file1, 'w')      # 0 lines
+        fd.close
+        fd = open(test_file2, 'w')      # 5 lines, all '\n'
+        for i in range(5):
+            fd.write('\n')
+        fd.close()
+        fd = open(test_file3, 'w')      # 25 chars, no \n, 1 lines
+        fd.write('no newline at end of line')
+        fd.close()
+        fd = open(test_file4, 'w')      # 1000 lines
+        for i in range(1000):
+            fd.write('The quick brown fox jumps over the lazy dog.\n')
+        fd.close()
+
+        # use file_length() to get and check lengths
+        size1 = file_length(test_file1)
+        msg = 'Expected file_length() to return 0, but got %d' % size1
+        self.failUnless(size1 == 0, msg)
+        size2 = file_length(test_file2)
+        msg = 'Expected file_length() to return 5, but got %d' % size2
+        self.failUnless(size2 == 5, msg)
+        size3 = file_length(test_file3)
+        msg = 'Expected file_length() to return 1, but got %d' % size3
+        self.failUnless(size3 == 1, msg)
+        size4 = file_length(test_file4)
+        msg = 'Expected file_length() to return 1000, but got %d' % size4
+        self.failUnless(size4 == 1000, msg)
+
 ################################################################################
 

branches/numpy/anuga/utilities/util_ext.h

@@ -346 +346 @@
 
 
-// check that numpy array objects arrays are C contiguous memory
+// check that numpy array objects are C contiguous memory
 #define CHECK_C_CONTIG(varname) if (!PyArray_ISCONTIGUOUS(varname)) { \
                                     char msg[1024]; \