Changeset 8578

Timestamp:

Sep 14, 2012, 9:56:39 PM (13 years ago)

Author:

steve

Message:

Added extra unit tests for set_stage and set_elevation operators

Location:

trunk/anuga_core/source/anuga

Files:

• __init__.py (modified) (2 diffs)
• abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (3 diffs)
• abstract_2d_finite_volumes/generic_domain.py (modified) (1 diff)
• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (1 diff)
• fit_interpolate/fit.py (modified) (21 diffs)
• fit_interpolate/general_fit_interpolate.py (modified) (2 diffs)
• operators/test_set_stage_operators.py (moved) (moved from trunk/anuga_core/source/anuga/operators/test_set_stage_operator.py) (2 diffs)
• pmesh/mesh_quadtree.py (modified) (5 diffs)
• shallow_water/boundaries.py (modified) (2 diffs)
• shallow_water/shallow_water_domain.py (modified) (14 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (1 diff)

trunk/anuga_core/source/anuga/__init__.py

@@ -73 +73 @@
 from anuga.shallow_water.boundaries import \
                     Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import \
+                    Compute_fluxes_boundary
 
 
@@ -172 +174 @@
     """
 
-    if 'verbose' in kwargs:
-        verbose = kwargs['verbose']
-        kwargs.pop('verbose')
-    else:
+    try:
+        verbose = kwargs.pop('verbose')
+    except:
         verbose = False
+
 
     points, vertices, boundary = rectangular_cross(*args, **kwargs)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -208 +208 @@
             Q.boundary_values[ids] = q_bdry[j]
 
+class Compute_fluxes_boundary(Boundary):
+    """ Associate the Compute_fluxes_boundary BC
+    to each boundary segment where a special calculation
+    of the fluxes will occur.
+    """
+
+    def __init__(self):
+        Boundary.__init__(self)
+        """ Instantiate a Compute_fluxes_boundary.
+            domain: underlying domain
+            """
+
+
+    def __repr__(self):
+        """ Return a representation of this instance. """
+        return 'Compute_fluxes_boundary(%s)' % self.domain
+
+    def evaluate(self, vol_id, edge_id):
+        """Do something basic"""
+
+        return None
+
+    def evaluate_segment(self, domain, segment_edges):
+        """Don't do any calculation when applying this BC.
+        The work will be done in compute_fluxes.
+        """
+
+        return
 
 
@@ -601 +629 @@
     def evaluate(self, vol_id=None, edge_id=None):
         """Return linearly interpolated values based on domain.time
-        at midpoint of segment defined by vol_id and edge_id.
+        at midpoint of segment defined by vol_id and edge_id.
         """
 
@@ -828 +856 @@
 
 
+

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

@@ -809 +809 @@
                 raise Exception(msg)
 
+
+
+        # Add a flag which can be used to distinguish flux boundaries within
+        # compute_fluxes_central
+        # Initialise to zero (which means 'not a flux_boundary')
+        self.boundary_flux_type = self.boundary_edges*0
+
+        # HACK to set the values of domain.boundary_flux
+        for k, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
+            # If Boundary set to Compute_fluxes_boundary identify as flux boundary
+            #print vol_id, edge_id, B
+
+            import anuga
+            if ( isinstance(B,anuga.Compute_fluxes_boundary) ):
+                self.boundary_flux_type[k]=1
+
+
     ##
     # @brief Set quantities based on a regional tag.

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -553 +553 @@
 
 
-        for j in range(self.boundary_length):
+        for j in xrange(self.boundary_length):
             id  = self.boundary_cells[j]
             edge = self.boundary_edges[j]

trunk/anuga_core/source/anuga/fit_interpolate/fit.py

@@ -46 +46 @@
 
 import numpy as num
+import sys
 
 
@@ -100 +101 @@
                  triangles,
                  mesh,
-                 mesh_origin,
-                 verbose)
+                 mesh_origin=mesh_origin,
+                 verbose=verbose)
         
         m = self.mesh.number_of_nodes # Nbr of basis functions (vertices)
@@ -110 +111 @@
         self.point_count = 0
         if self.alpha <> 0:
-            if verbose: log.critical('Building smoothing matrix')
-            self._build_smoothing_matrix_D()
-            
+            if verbose: log.critical('Fit: Building smoothing matrix')
+            self._build_smoothing_matrix_D(verbose=verbose)
+
+        if verbose: log.critical('Fit: Get Boundary Polygon')
         bd_poly = self.mesh.get_boundary_polygon()    
         self.mesh_boundary_polygon = ensure_numeric(bd_poly)
+
+        if verbose: log.critical('Fit: Done')
+
             
     def _build_coefficient_matrix_B(self,
@@ -126 +131 @@
         """
 
+        if verbose:
+            print 'Fit: Build Coefficient Matrix B'
+
+
         if self.alpha <> 0:
             #if verbose: log.critical('Building smoothing matrix')
             #self._build_smoothing_matrix_D()
+            # FIXME SR: This is quite time consuming.
+            # As AtA and D have same structure it should be possible
+            # to speed this up.
             self.B = self.AtA + self.alpha*self.D
         else:
             self.B = self.AtA
 
+
+        if verbose:
+            print 'Fit: Convert Coefficient Matrix B to Sparse_CSR'
+            
         # Convert self.B matrix to CSR format for faster matrix vector
         self.B = Sparse_CSR(self.B)
 
-    def _build_smoothing_matrix_D(self):
+    def _build_smoothing_matrix_D(self, verbose=False):
         """Build m x m smoothing matrix, where
         m is the number of basis functions phi_k (one per vertex)
@@ -169 +185 @@
         self.D = Sparse(m,m)
 
+        if verbose :
+            print '['+60*' '+']',
+            sys.stdout.flush()
+
+        N = len(self.mesh)
+        M = N/60
+
         # For each triangle compute contributions to D = D1+D2
-        for i in range(len(self.mesh)):
+        for i in xrange(N):
+
+            if verbose and i%M == 0 :
+                #restart_line()
+                print '\r['+(i/M)*'.'+(60-(i/M))*' ' +']',
+                sys.stdout.flush()
 
             # Get area
@@ -212 +240 @@
             self.D[v2,v0] += e20
             self.D[v0,v2] += e20
+
+        if verbose:
+            print ''
 
     def get_D(self):
@@ -242 +273 @@
         The number of attributes of the data points does not change
         """
-        
+
+
+        if verbose:
+            print 'Fit: Build Matrix AtA Atz'
+
         # Build n x m interpolation matrix
         if self.AtA == None:
@@ -271 +306 @@
         # Compute matrix elements for points inside the mesh
         triangles = self.mesh.triangles # Shorthand
-        for d, i in enumerate(inside_indices):
+
+
+        if verbose :
+            print '['+60*' '+']',
+            sys.stdout.flush()
+
+        m = n/60
+
+
+        for i in xrange(n):
+            d = inside_indices[i]
             # For each data_coordinate point
             # if verbose and d%((n+10)/10)==0: log.critical('Doing %d of %d'
                                                             # %(d, n))
-            x = point_coordinates[i]
+
+
+            if verbose and i%m == 0 :
+                print '\r['+(i/m)*'.'+(60-(i/m))*' '+']',
+                sys.stdout.flush()
+
+
+            x = point_coordinates[d]
             
             element_found, sigma0, sigma1, sigma2, k = \
@@ -302 +354 @@
                 
                 # data point has fallen within a hole - so ignore it.
-                
+
+
+        if verbose:
+            print ''
+            
         self.AtA = AtA
 
@@ -324 +380 @@
           
         """
-        
+
+
+        if verbose:
+            print 'Fit.fit: Initializing'
+
         # Use blocking to load in the point info
         if isinstance(point_coordinates_or_filename, basestring):
@@ -357 +417 @@
                 points = geo_block.get_data_points(absolute=True)
                 z = geo_block.get_attributes(attribute_name=attribute_name)
-                self.build_fit_subset(points, z, verbose=verbose)
+                self.build_fit_subset(points, z, attribute_name, verbose)
 
                 # FIXME(Ole): I thought this test would make sense here
@@ -368 +428 @@
         else:
             point_coordinates =  point_coordinates_or_filename
-            
+
+
+
         if point_coordinates is None:
-            if verbose: log.critical('Warning: no data points in fit')
+            if verbose: log.critical('Fit.fit: Warning: no data points in fit')
             msg = 'No interpolation matrix.'
             assert self.AtA is not None, msg
@@ -381 +443 @@
             # if isinstance(point_coordinates,Geospatial_data) and z is None:
             # z will come from the geo-ref
-            self.build_fit_subset(point_coordinates, z, verbose)
+            self.build_fit_subset(point_coordinates, z, verbose=verbose)
+
+
+
+
 
         # Check sanity
@@ -393 +459 @@
             msg += 'positive value,\ne.g. 1.0e-3.'
             raise TooFewPointsError(msg)
+
 
         self._build_coefficient_matrix_B(verbose)
@@ -408 +475 @@
             #raise VertsWithNoTrianglesError(msg)
         
-        
+        if verbose:
+            print 'Fit.fit: Solve Fitting Equation'
+
         return conjugate_gradient(self.B, self.Atz, self.Atz,
                                   imax=2*len(self.Atz) )
@@ -575 +644 @@
 
         if verbose: log.critical('FitInterpolate: Building mesh')
-        mesh = Mesh(vertex_coordinates, triangles)
-        mesh.check_integrity()
+        mesh = Mesh(vertex_coordinates, triangles, verbose=verbose)
+        #mesh.check_integrity() # expensive
     
     
@@ -648 +717 @@
         old_title_list = mesh_dict['vertex_attribute_titles']
 
-    if verbose: log.critical('tsh file %s loaded' % mesh_file)
+    if verbose: log.critical('Fit_to_Mesh_File: tsh file %s loaded' % mesh_file)
 
     # load in the points file
@@ -668 +737 @@
         mesh_origin = None
 
-    if verbose: log.critical("points file loaded")
-    if verbose: log.critical("fitting to mesh")
+    if verbose: log.critical("Fit_to_Mesh_File: points file loaded")
+    if verbose: log.critical("Fit_to_Mesh_File: fitting to mesh")
     f = fit_to_mesh(point_coordinates,
                     vertex_coordinates,
@@ -679 +748 @@
                     data_origin = None,
                     mesh_origin = mesh_origin)
-    if verbose: log.critical("finished fitting to mesh")
+    if verbose: log.critical("Fit_to_Mesh_File: finished fitting to mesh")
 
     # convert array to list of lists
@@ -688 +757 @@
         old_title_list.extend(title_list)
         #FIXME can this be done a faster way? - DSG
-        for i in range(len(old_point_attributes)):
+        for i in xrange(len(old_point_attributes)):
             old_point_attributes[i].extend(new_point_attributes[i])
         mesh_dict['vertex_attributes'] = old_point_attributes
@@ -696 +765 @@
         mesh_dict['vertex_attribute_titles'] = title_list
 
-    if verbose: log.critical("exporting to file %s" % mesh_output_file)
+    if verbose: log.critical("Fit_to_Mesh_File: exporting to file %s" % mesh_output_file)
 
     try:

trunk/anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

@@ -91 +91 @@
 
                 if verbose: log.critical('FitInterpolate: Building mesh')
-                self.mesh = Mesh(vertex_coordinates, triangles)
+                self.mesh = Mesh(vertex_coordinates, triangles, verbose=verbose)
                 #self.mesh.check_integrity() # Time consuming
             else:
@@ -102 +102 @@
             #This stores indices of vertices
             t0 = time.time()
-            self.root = MeshQuadtree(self.mesh)
+            self.root = MeshQuadtree(self.mesh, verbose=verbose)
         
             build_quadtree_time =  time.time()-t0

trunk/anuga_core/source/anuga/operators/test_set_stage_operators.py

@@ -21 +21 @@
 
 
-class Test_set_operators(unittest.TestCase):
+class Test_set_stage_operators(unittest.TestCase):
     def setUp(self):
         pass
@@ -142 +142 @@
 
 
-
-
-            
+    def test_set_stage_operator_function(self):
+        from anuga.config import rho_a, rho_w, eta_w
+        from math import pi, cos, sin
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+
+
+
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        # Apply operator to these triangles
+        indices = [0,1,3]
+
+
+        def stage(t):
+            if t < 10.0:
+                return 5.0
+            else:
+                return 10.0
+
+        operator = Set_stage_operator(domain, stage=stage, indices=indices)
+
+        # Apply Operator at time t=1.0
+        domain.set_time(1.0)
+        operator()
+
+        stage_ex = [ 5.,  5.,   1.,  5.]
+
+
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
+        # Apply Operator at time t=15.0
+        domain.set_time(15.0)
+        operator()
+
+        stage_ex = [ 10.,  10.,   1.,  10.]
+
+
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
+
+    def test_set_stage_operator_large_function(self):
+        from anuga.config import rho_a, rho_w, eta_w
+        from math import pi, cos, sin
+
+
+        length = 2.0
+        width = 2.0
+        dx = dy = 0.5
+        domain = rectangular_cross_domain(int(length/dx), int(width/dy),
+                                              len1=length, len2=width)
+
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0.0)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0)
+
+        R = Reflective_boundary(domain)
+        domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
+
+
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+
+
+        def stage(t):
+            if t < 10.0:
+                return 5.0
+            else:
+                return 7.0
+
+        polygon = [(0.5,0.5), (1.5,0.5), (1.5,1.5), (0.5,1.5)]
+        operator = Polygonal_set_stage_operator(domain, stage=stage, polygon=polygon)
+
+        # Apply Operator at time t=1.0
+        domain.set_time(1.0)
+        operator()
+
+        stage_ex = [ 1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+                     1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+                     5.0,  5.0,  5.0,  5.0,  5.0,  5.0,  5.0,  5.0,  1.0,  1.0,
+                     1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  5.0,  5.0,  5.0,  5.0,
+                     5.0,  5.0,  5.0,  5.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+                     1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+                     1.0,  1.0,  1.0,  1.0]
+
+
+#        print domain.quantities['elevation'].centroid_values
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
+        # Apply Operator at time t=15.0
+        domain.set_time(15.0)
+        operator()
+
+        stage_ex = [ 1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+                     1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+                     7.0,  7.0,  7.0,  7.0,  7.0,  7.0,  7.0,  7.0,  1.0,  1.0,
+                     1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  7.0,  7.0,  7.0,  7.0,
+                     7.0,  7.0,  7.0,  7.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+                     1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+                     1.0,  1.0,  1.0,  1.0]
+
+
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
+
+
+
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_set_operators, 'test')
+    suite = unittest.makeSuite(Test_set_stage_operators, 'test')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)

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

@@ -24 +24 @@
 
 import numpy as num
+import sys
  
 
@@ -41 +42 @@
         It contains optimisations and search patterns specific to meshes.
     """
-    def __init__(self, mesh):
+    def __init__(self, mesh, verbose=False):
         """Build quad tree for mesh.
 
@@ -60 +61 @@
         
         normals = mesh.get_normals()
-        
+
+        if verbose :
+            print '['+60*' '+']',
+            sys.stdout.flush()
+
+        M = N/60
+
         # Check each triangle
-        for i in range(N):
+        for i in xrange(N):
+
+            if verbose and i%M == 0 :
+                #restart_line()
+                print '\r['+(i/M)*'.'+(60-(i/M))*' ' +']',
+                sys.stdout.flush()
+
             i3 = 3*i
             x0, y0 = V[i3, :]
@@ -68 +81 @@
             x2, y2 = V[i3+2, :]
 
+            #FIXME SR: Should save memory by just using triangle id!
             node_data = [i, V[i3:i3+3, :], normals[i, :]]
 
@@ -74 +88 @@
                              min([y0, y1, y2]), max([y0, y1, y2])), \
                              node_data))
+
+        if verbose:
+            print ''
 
     def search_fast(self, point):

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

@@ -477 +477 @@
             
         # Assign conserved quantities and return
-        q = num.array([elevation + depth, xmomentum, ymomentum], num.Float)
+        q = num.array([elevation + depth, xmomentum, ymomentum], num.float)
         return q
 
@@ -580 +580 @@
                 q[j] += self.mean_stage
         return q
+
+
+
+
+

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

@@ -1042 +1042 @@
 
         elif self.compute_fluxes_method == 'tsunami':
-            # Using Gareth Davies well nbalanced scheme
+            # Using Gareth Davies well balanced scheme
             # Flux calculation and gravity incorporated in same
             # procedure
@@ -1054 +1054 @@
 
             # Shortcuts
-            Stage = domain.quantities['stage']
-            Xmom = domain.quantities['xmomentum']
-            Ymom = domain.quantities['ymomentum']
-            Bed = domain.quantities['elevation']
+            Stage = self.quantities['stage']
+            Xmom = self.quantities['xmomentum']
+            Ymom = self.quantities['ymomentum']
+            Bed = self.quantities['elevation']
 
             timestep = self.evolve_max_timestep
@@ -1109 +1109 @@
         if self.compute_fluxes_method == 'tsunami':
 
+
+            # FIXME SR: Clean up code to just take self (domain) as
+            # input argument
             from swb2_domain_ext import protect
 
@@ -1120 +1123 @@
             areas = self.areas
 
-            protect(self.minimum_allowed_height, domain.maximum_allowed_speed,
-                domain.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas)
-
-
-            from swb2_domain_ext import extrapolate_second_order_edge_sw as extrapol2
+            protect(self.minimum_allowed_height, self.maximum_allowed_speed,
+                self.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas)
+
+
+            from swb2_domain_ext import extrapolate_second_order_edge_sw as extrapol2_ext
 
             # Shortcuts
@@ -1132 +1135 @@
             Elevation = self.quantities['elevation']
 
-            extrapol2(self,
+            extrapol2_ext(self,
                   self.surrogate_neighbours,
                   self.number_of_boundaries,
@@ -1154 +1157 @@
 
         else:
+            # Code for original method
             if self.use_edge_limiter:
                 self.distribute_using_edge_limiter()
             else:
-                #distribute_using_vertex_limiter(self)
                 self.distribute_using_vertex_limiter()
 
@@ -1193 +1196 @@
                 raise Exception('Unknown order')
 
-        balance_deep_and_shallow(self)
+        self.balance_deep_and_shallow()
 
         # Compute edge values by interpolation
@@ -1200 +1203 @@
             Q.interpolate_from_vertices_to_edges()
 
-    def distribute_using_vertex_limiter(domain):
+    def distribute_using_vertex_limiter(self):
         """Distribution from centroids to vertices specific to the SWW equation.
 
@@ -1221 +1224 @@
 
         # Remove very thin layers of water
-        domain.protect_against_infinitesimal_and_negative_heights()
+        self.protect_against_infinitesimal_and_negative_heights()
 
         # Extrapolate all conserved quantities
-        if domain.optimised_gradient_limiter:
+        if self.optimised_gradient_limiter:
             # MH090605 if second order,
             # perform the extrapolation and limiting on
             # all of the conserved quantities
 
-            if (domain._order_ == 1):
-                for name in domain.conserved_quantities:
-                    Q = domain.quantities[name]
+            if (self._order_ == 1):
+                for name in self.conserved_quantities:
+                    Q = self.quantities[name]
                     Q.extrapolate_first_order()
-            elif domain._order_ == 2:
-                domain.extrapolate_second_order_sw()
+            elif self._order_ == 2:
+                self.extrapolate_second_order_sw()
             else:
                 raise Exception('Unknown order')
@@ -1240 +1243 @@
             # Old code:
             for name in domain.conserved_quantities:
-                Q = domain.quantities[name]
-
-                if domain._order_ == 1:
+                Q = self.quantities[name]
+
+                if self._order_ == 1:
                     Q.extrapolate_first_order()
-                elif domain._order_ == 2:
+                elif self._order_ == 2:
                     Q.extrapolate_second_order_and_limit_by_vertex()
                 else:
@@ -1250 +1253 @@
 
         # Take bed elevation into account when water heights are small
-        balance_deep_and_shallow(domain)
+        self.balance_deep_and_shallow()
 
         # Compute edge values by interpolation
-        for name in domain.conserved_quantities:
-            Q = domain.quantities[name]
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
             Q.interpolate_from_vertices_to_edges()
 
@@ -1267 +1270 @@
 
             # shortcuts
-            wc = domain.quantities['stage'].centroid_values
-            wv = domain.quantities['stage'].vertex_values
-            zc = domain.quantities['elevation'].centroid_values
-            zv = domain.quantities['elevation'].vertex_values
-            xmomc = domain.quantities['xmomentum'].centroid_values
-            ymomc = domain.quantities['ymomentum'].centroid_values
-            areas = domain.areas
-
-            protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
-                domain.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas)
+            wc = self.quantities['stage'].centroid_values
+            wv = self.quantities['stage'].vertex_values
+            zc = self.quantities['elevation'].centroid_values
+            zv = self.quantities['elevation'].vertex_values
+            xmomc = self.quantities['xmomentum'].centroid_values
+            ymomc = self.quantities['ymomentum'].centroid_values
+            areas = self.areas
+
+            protect(self.minimum_allowed_height, self.maximum_allowed_speed,
+                self.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas)
 
         else:
@@ -1290 +1293 @@
                     self.epsilon, wc, zc, xmomc, ymomc)
 
-        
+
+    def balance_deep_and_shallow(self):
+        """Compute linear combination between stage as computed by
+        gradient-limiters limiting using w, and stage computed by
+        gradient-limiters limiting using h (h-limiter).
+        The former takes precedence when heights are large compared to the
+        bed slope while the latter takes precedence when heights are
+        relatively small.  Anything in between is computed as a balanced
+        linear combination in order to avoid numerical disturbances which
+        would otherwise appear as a result of hard switching between
+        modes.
+
+        Wrapper for C implementation
+        """
+
+        from shallow_water_ext import balance_deep_and_shallow \
+                                      as balance_deep_and_shallow_ext
+
+        # Shortcuts
+        wc = self.quantities['stage'].centroid_values
+        zc = self.quantities['elevation'].centroid_values
+        wv = self.quantities['stage'].vertex_values
+        zv = self.quantities['elevation'].vertex_values
+
+        # Momentums at centroids
+        xmomc = self.quantities['xmomentum'].centroid_values
+        ymomc = self.quantities['ymomentum'].centroid_values
+
+        # Momentums at vertices
+        xmomv = self.quantities['xmomentum'].vertex_values
+        ymomv = self.quantities['ymomentum'].vertex_values
+
+        balance_deep_and_shallow_ext(self,
+                                   wc, zc, wv, zv, wc,
+                                   xmomc, ymomc, xmomv, ymomv)
 
     def update_other_quantities(self):
@@ -1990 +2027 @@
 ##            domain.epsilon, wc, zc, xmomc, ymomc)
 
-def balance_deep_and_shallow(domain):
-    """Compute linear combination between stage as computed by
-    gradient-limiters limiting using w, and stage computed by
-    gradient-limiters limiting using h (h-limiter).
-    The former takes precedence when heights are large compared to the
-    bed slope while the latter takes precedence when heights are
-    relatively small.  Anything in between is computed as a balanced
-    linear combination in order to avoid numerical disturbances which
-    would otherwise appear as a result of hard switching between
-    modes.
-
-    Wrapper for C implementation
-    """
-
-    from anuga.shallow_water.shallow_water_ext import balance_deep_and_shallow \
-                                  as balance_deep_and_shallow_c
-
-    # Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    wv = domain.quantities['stage'].vertex_values
-    zv = domain.quantities['elevation'].vertex_values
-
-    # Momentums at centroids
-    xmomc = domain.quantities['xmomentum'].centroid_values
-    ymomc = domain.quantities['ymomentum'].centroid_values
-
-    # Momentums at vertices
-    xmomv = domain.quantities['xmomentum'].vertex_values
-    ymomv = domain.quantities['ymomentum'].vertex_values
-
-    balance_deep_and_shallow_c(domain,
-                               wc, zc, wv, zv, wc,
-                               xmomc, ymomc, xmomv, ymomv)
+
 
 

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

@@ -7012 +7012 @@
         # Large test set revealed one problem
         points_file = os.path.join(path, 'test_points_large.csv')
+
+        domain.set_quantity('elevation', filename=points_file,
+                                use_cache=False, verbose=verbose)
         try:
             domain.set_quantity('elevation', filename=points_file,