Changeset 5892

Timestamp:

Nov 5, 2008, 4:29:38 PM (16 years ago)

Author:

rwilson

Message:

Made numpy changes, lots of problems.

Location:

anuga_core/source/anuga/abstract_2d_finite_volumes

Files:

• domain.py (modified) (15 diffs)
• ermapper_grids.py (modified) (9 diffs)
• general_mesh.py (modified) (17 diffs)
• generic_boundary_conditions.py (modified) (7 diffs)
• mesh_factory.py (modified) (9 diffs)
• neighbour_mesh.py (modified) (20 diffs)
• pmesh2domain.py (modified) (3 diffs)
• quantity.py (modified) (24 diffs)
• quantity_ext.c (modified) (1 diff)
• region.py (modified) (2 diffs)
• show_balanced_limiters.py (modified) (2 diffs)
• test_domain.py (modified) (18 diffs)
• test_ermapper.py (modified) (8 diffs)
• test_general_mesh.py (modified) (20 diffs)
• test_generic_boundary_conditions.py (modified) (5 diffs)
• test_ghost.py (modified) (2 diffs)
• test_neighbour_mesh.py (modified) (51 diffs)
• test_pmesh2domain.py (modified) (6 diffs)
• test_quantity.py (modified) (124 diffs)
• test_region.py (modified) (14 diffs)
• test_util.py (modified) (67 diffs)
• util.py (modified) (20 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -8 +8 @@
 """
 
-from Numeric import allclose, argmax, zeros, Float
+##from numpy.oldnumeric import allclose, argmax, zeros, Float
+import numpy
+
 from anuga.config import epsilon
 from anuga.config import beta_euler, beta_rk2
@@ -105 +107 @@
 
         if verbose: print 'Initialising Domain'
-        from Numeric import zeros, Float, Int, ones
         from quantity import Quantity
 
@@ -155 +156 @@
         for key in self.full_send_dict:
             buffer_shape = self.full_send_dict[key][0].shape[0]
-            self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
+            self.full_send_dict[key].append(numpy.zeros( (buffer_shape,self.nsys), numpy.float))
 
 
         for key in self.ghost_recv_dict:
             buffer_shape = self.ghost_recv_dict[key][0].shape[0]
-            self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
+            self.ghost_recv_dict[key].append(numpy.zeros( (buffer_shape,self.nsys), numpy.float))
 
 
@@ -168 +169 @@
         N = len(self) #number_of_elements
         self.number_of_elements = N
-        self.tri_full_flag = ones(N, Int)
+        self.tri_full_flag = numpy.ones(N, numpy.int)
         for i in self.ghost_recv_dict.keys():
             for id in self.ghost_recv_dict[i][0]:
@@ -175 +176 @@
         # Test the assumption that all full triangles are store before
         # the ghost triangles.
-        if not allclose(self.tri_full_flag[:self.number_of_full_nodes],1):
+        if not numpy.allclose(self.tri_full_flag[:self.number_of_full_nodes],1):
             print 'WARNING:  Not all full triangles are store before ghost triangles'
                         
@@ -230 +231 @@
         # calculation
         N = len(self) # Number_of_triangles
-        self.already_computed_flux = zeros((N, 3), Int)
+        self.already_computed_flux = numpy.zeros((N, 3), numpy.int)
 
         # Storage for maximal speeds computed for each triangle by
         # compute_fluxes
         # This is used for diagnostics only (reset at every yieldstep)
-        self.max_speed = zeros(N, Float)
+        self.max_speed = numpy.zeros(N, numpy.float)
 
         if mesh_filename is not None:
@@ -266 +267 @@
         """
 
-        from Numeric import zeros, Float
-
         if not (vertex is None or edge is None):
             msg = 'Values for both vertex and edge was specified.'
@@ -273 +272 @@
             raise msg
 
-        q = zeros( len(self.conserved_quantities), Float)
+        q = numpy.zeros( len(self.conserved_quantities), numpy.float)
 
         for i, name in enumerate(self.conserved_quantities):
@@ -769 +768 @@
             # Find index of largest computed flux speed
             if triangle_id is None:
-                k = self.k = argmax(self.max_speed)
+                k = self.k = numpy.argmax(self.max_speed)
             else:
                 errmsg = 'Triangle_id %d does not exist in mesh: %s' %(triangle_id,
@@ -1086 +1085 @@
 
         """
+        print '.evolve: 0'
 
         from anuga.config import min_timestep, max_timestep, epsilon
@@ -1096 +1096 @@
                %self.get_boundary_tags()
         assert hasattr(self, 'boundary_objects'), msg
+        print '.evolve: 1'
 
 
@@ -1104 +1105 @@
 
         self._order_ = self.default_order
+        print '.evolve: 2'
 
 
@@ -1127 +1129 @@
         self.number_of_first_order_steps = 0
 
-
+        print '.evolve: 3'
         # Update ghosts
         self.update_ghosts()
 
+        print '.evolve: 4'
         # Initial update of vertex and edge values
         self.distribute_to_vertices_and_edges()
 
+        print '.evolve: 5'
         # Update extrema if necessary (for reporting)
         self.update_extrema()
         
+        print '.evolve: 6'
         # Initial update boundary values
         self.update_boundary()
 
+        print '.evolve: 7'
         # Or maybe restore from latest checkpoint
         if self.checkpoint is True:
             self.goto_latest_checkpoint()
+        print '.evolve: 8'
 
         if skip_initial_step is False:
@@ -1202 +1209 @@
                 self.number_of_steps = 0
                 self.number_of_first_order_steps = 0
-                self.max_speed = zeros(N, Float)
+                self.max_speed = numpy.zeros(N, numpy.float)
 
     def evolve_one_euler_step(self, yieldstep, finaltime):
@@ -1568 +1575 @@
         """
 
-        from Numeric import ones, sum, equal, Float
-
         N = len(self) # Number_of_triangles
         d = len(self.conserved_quantities)

anuga_core/source/anuga/abstract_2d_finite_volumes/ermapper_grids.py

@@ -2 +2 @@
 
 # from os import open, write, read
-import Numeric
-
-celltype_map = {'IEEE4ByteReal': Numeric.Float32, 'IEEE8ByteReal': Numeric.Float64}
+##import numpy.oldnumeric as Numeric
+import numpy
+
+celltype_map = {'IEEE4ByteReal': numpy.float32, 'IEEE8ByteReal': numpy.float64}
 
 
@@ -11 +12 @@
     write_ermapper_grid(ofile, data, header = {}):
     
-    Function to write a 2D Numeric array to an ERMapper grid.  There are a series of conventions adopted within
+    Function to write a 2D NumPy array to an ERMapper grid.  There are a series of conventions adopted within
     this code, specifically:
     1)  The registration coordinate for the data is the SW (or lower-left) corner of the data
     2)  The registration coordinates refer to cell centres
-    3)  The data is a 2D Numeric array with the NW-most data in element (0,0) and the SE-most data in element (N,M)
+    3)  The data is a 2D NumPy array with the NW-most data in element (0,0) and the SE-most data in element (N,M)
         where N is the last line and M is the last column
     4)  There has been no testng of the use of a rotated grid.  Best to keep data in an NS orientation
@@ -22 +23 @@
     ofile:      string - filename for output (note the output will consist of two files
                 ofile and ofile.ers.  Either of these can be entered into this function
-    data:       Numeric.array - 2D array containing the data to be output to the grid
+    data:       NumPy.array - 2D array containing the data to be output to the grid
     header:     dictionary - contains spatial information about the grid, in particular:
                     header['datum'] datum for the data ('"GDA94"')
@@ -57 +58 @@
 
     # Check that the data is a 2 dimensional array
-    data_size = Numeric.shape(data)
+    data_size = numpy.shape(data)
     assert len(data_size) == 2
     
@@ -83 +84 @@
     nrofcellsperlines = int(header['nrofcellsperline'])
     data = read_ermapper_data(data_file)
-    data = Numeric.reshape(data,(nroflines,nrofcellsperlines))
+    data = numpy.reshape(data,(nroflines,nrofcellsperlines))
     return data
     
@@ -163 +164 @@
     return header                      
 
-def write_ermapper_data(grid, ofile, data_format = Numeric.Float32):
+def write_ermapper_data(grid, ofile, data_format = numpy.float32):
 
 
@@ -181 +182 @@
         
     
-    # Convert the array to data_format (default format is Float32)
+    # Convert the array to data_format (default format is float32)
     grid_as_float = grid.astype(data_format)
 
@@ -193 +194 @@
 
 
-def read_ermapper_data(ifile, data_format = Numeric.Float32):
+def read_ermapper_data(ifile, data_format = numpy.float32):
     # open input file in a binary format and read the input string
     fid = open(ifile,'rb')
@@ -199 +200 @@
     fid.close()
 
-    # convert input string to required format (Note default format is Numeric.Float32)
-    grid_as_float = Numeric.fromstring(input_string,data_format)
+    # convert input string to required format (Note default format is numpy.float32)
+    grid_as_float = numpy.fromstring(input_string,data_format)
     return grid_as_float
 

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -1 @@
-
-from Numeric import concatenate, reshape, take, allclose
-from Numeric import array, zeros, Int, Float, sqrt, sum, arange
+import numpy
 
 from anuga.coordinate_transforms.geo_reference import Geo_reference
@@ -90 +88 @@
         if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' 
 
-        self.triangles = array(triangles, Int)
-        self.nodes = array(nodes, Float)
+##        self.triangles = array(triangles, Int)
+##        self.nodes = array(nodes, Float)
+        self.triangles = numpy.array(triangles, numpy.int)
+        self.nodes = numpy.array(nodes, numpy.float)
 
 
@@ -138 +138 @@
 
         msg = 'Vertex indices reference non-existing coordinate sets'
-        assert max(self.triangles.flat) < self.nodes.shape[0], msg
+        assert max(self.triangles.ravel()) < self.nodes.shape[0], msg
 
 
@@ -146 +146 @@
                       max(self.nodes[:,0]), max(self.nodes[:,1]) ]
 
-        self.xy_extent = array(xy_extent, Float)
+##        self.xy_extent = array(xy_extent, Float)
+        self.xy_extent = numpy.array(xy_extent, numpy.float)
 
 
         # Allocate space for geometric quantities
-        self.normals = zeros((N, 6), Float)
-        self.areas = zeros(N, Float)
-        self.edgelengths = zeros((N, 3), Float)
+##        self.normals = zeros((N, 6), Float)
+##        self.areas = zeros(N, Float)
+##        self.edgelengths = zeros((N, 3), Float)
+        self.normals = numpy.zeros((N, 6), numpy.float)
+        self.areas = numpy.zeros(N, numpy.float)
+        self.edgelengths = numpy.zeros((N, 3), numpy.float)
 
         # Get x,y coordinates for all triangles and store
@@ -187 +191 @@
             #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
 
-            n0 = array([x2 - x1, y2 - y1])
-            l0 = sqrt(sum(n0**2))
-
-            n1 = array([x0 - x2, y0 - y2])
-            l1 = sqrt(sum(n1**2))
-
-            n2 = array([x1 - x0, y1 - y0])
-            l2 = sqrt(sum(n2**2))
+            n0 = numpy.array([x2 - x1, y2 - y1])
+            l0 = numpy.sqrt(numpy.sum(n0**2))
+
+            n1 = numpy.array([x0 - x2, y0 - y2])
+            l1 = numpy.sqrt(numpy.sum(n1**2))
+
+            n2 = numpy.array([x1 - x0, y1 - y0])
+            l2 = numpy.sqrt(numpy.sum(n2**2))
 
             # Normalise
@@ -274 +278 @@
         if absolute is True:
             if not self.geo_reference.is_absolute():
-                return V + array([self.geo_reference.get_xllcorner(),
-                                  self.geo_reference.get_yllcorner()])
+                return V + numpy.array([self.geo_reference.get_xllcorner(),
+                                        self.geo_reference.get_yllcorner()])
             else:
                 return V
@@ -311 +315 @@
             i = triangle_id
             msg = 'triangle_id must be an integer'
+            print 'type(triangle_id)=%s. triangle_id=%s' % (type(triangle_id), str(triangle_id))
             assert int(i) == i, msg
             assert 0 <= i < self.number_of_triangles
@@ -316 +321 @@
             i3 = 3*i  
             if absolute is True and not self.geo_reference.is_absolute():
-                offset=array([self.geo_reference.get_xllcorner(),
-                                  self.geo_reference.get_yllcorner()])
-                return array([V[i3,:]+offset,
-                              V[i3+1,:]+offset,
-                              V[i3+2,:]+offset])
+                offset = numpy.array([self.geo_reference.get_xllcorner(),
+                                      self.geo_reference.get_yllcorner()])
+                return numpy.array([V[i3,:]+offset,
+                                    V[i3+1,:]+offset,
+                                    V[i3+2,:]+offset])
             else:
-                return array([V[i3,:], V[i3+1,:], V[i3+2,:]])
+                return numpy.array([V[i3,:], V[i3+1,:], V[i3+2,:]])
                 
 
@@ -349 +354 @@
 
         M = self.number_of_triangles
-        vertex_coordinates = zeros((3*M, 2), Float)
+##        vertex_coordinates = zeros((3*M, 2), Float)
+        vertex_coordinates = numpy.zeros((3*M, 2), numpy.float)
 
         for i in range(M):
@@ -376 +382 @@
             indices = range(M)
 
-        return take(self.triangles, indices)
+        return numpy.take(self.triangles, indices, axis=0)
     
 
@@ -409 +415 @@
         
         #T = reshape(array(range(K)).astype(Int), (M,3))
-        T = reshape(arange(K).astype(Int), (M,3))  # Faster
+##        T = reshape(arange(K).astype(Int), (M,3))  # Faster
+        T = numpy.reshape(numpy.arange(K).astype(numpy.int), (M,3))  # Faster
         
         return T     
@@ -439 +446 @@
         if node is not None:
             # Get index for this node
-            first = sum(self.number_of_triangles_per_node[:node])
+            first = numpy.sum(self.number_of_triangles_per_node[:node])
             
             # Get number of triangles for this node
-            count = self.number_of_triangles_per_node[node]
+            count = int(self.number_of_triangles_per_node[node])
 
             for i in range(count):
@@ -452 +459 @@
                 triangle_list.append( (volume_id, vertex_id) )
 
-            triangle_list = array(triangle_list)    
+            triangle_list = numpy.array(triangle_list)    
         else:
             # Get info for all nodes recursively.
@@ -529 +536 @@
 
         # Count number of triangles per node
-        number_of_triangles_per_node = zeros(self.number_of_full_nodes)
+        number_of_triangles_per_node = numpy.zeros(self.number_of_full_nodes)
         for volume_id, triangle in enumerate(self.get_triangles()):
             for vertex_id in triangle:
@@ -535 +542 @@
 
         # Allocate space for inverted structure
-        number_of_entries = sum(number_of_triangles_per_node)
-        vertex_value_indices = zeros(number_of_entries)
+        number_of_entries = numpy.sum(number_of_triangles_per_node)
+        vertex_value_indices = numpy.zeros(number_of_entries)
 
         # Register (triangle, vertex) indices for each node
@@ -582 +589 @@
         Y = C[:,1:6:2].copy()
 
-        xmin = min(X.flat)
-        xmax = max(X.flat)
-        ymin = min(Y.flat)
-        ymax = max(Y.flat)
+        xmin = min(X.ravel())
+        xmax = max(X.ravel())
+        ymin = min(Y.ravel())
+        ymax = max(Y.ravel())
         #print "C",C
         return xmin, xmax, ymin, ymax
@@ -598 +605 @@
         """
 
-        return sum(self.areas)
-
-        
-        
+        return numpy.sum(self.areas)

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -7 +7 @@
 from anuga.fit_interpolate.interpolate import Modeltime_too_late
 from anuga.fit_interpolate.interpolate import Modeltime_too_early
+
+import numpy
 
 
@@ -67 +69 @@
             raise Exception, msg
 
-        from Numeric import array, Float
-        self.conserved_quantities=array(conserved_quantities).astype(Float)
+        self.conserved_quantities=numpy.array(conserved_quantities).astype(numpy.float)
 
     def __repr__(self):
@@ -96 +97 @@
             raise msg
 
-
-        from Numeric import array, Float
         try:
-            q = array(q).astype(Float)
+            q = numpy.array(q).astype(numpy.float)
         except:
             msg = 'Return value from time boundary function could '
@@ -136 +135 @@
     def __init__(self, filename, domain):
         import time
-        from Numeric import array
         from anuga.config import time_format
         from anuga.abstract_2d_finite_volumes.util import File_function
@@ -206 +204 @@
 
         import time
-        from Numeric import array, zeros, Float
         from anuga.config import time_format
         from anuga.abstract_2d_finite_volumes.util import file_function
@@ -222 +219 @@
 
         if verbose: print 'Find midpoint coordinates of entire boundary'
-        self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
+        self.midpoint_coordinates = numpy.zeros( (len(domain.boundary), 2), numpy.float)
         boundary_keys = domain.boundary.keys()
 
@@ -242 +239 @@
             
             # Compute midpoints
-            if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
-            if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
-            if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])
+            if edge_id == 0: m = numpy.array([(x1 + x2)/2, (y1 + y2)/2])
+            if edge_id == 1: m = numpy.array([(x0 + x2)/2, (y0 + y2)/2])
+            if edge_id == 2: m = numpy.array([(x1 + x0)/2, (y1 + y0)/2])
 
             # Convert to absolute UTM coordinates

anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py

@@ -2 +2 @@
 mesh file formats
 """
+
+import numpy
 
 
@@ -73 +75 @@
 
     from anuga.config import epsilon
-    from Numeric import zeros, Float, Int
 
     delta1 = float(len1)/m
@@ -93 +94 @@
     index = Index(n,m)
 
-    points = zeros( (Np,2), Float)
+    points = numpy.zeros( (Np,2), numpy.float)
 
     for i in range(m+1):
@@ -105 +106 @@
 
 
-    elements = zeros( (Nt,3), Int)
+    elements = numpy.zeros( (Nt,3), numpy.int)
     boundary = {}
     nt = -1
@@ -149 +150 @@
 
     from anuga.config import epsilon
-    from Numeric import zeros, Float, Int
 
     delta1 = float(len1)/m
@@ -208 +208 @@
     """
 
-    from Numeric import array
     import math
 
@@ -510 +509 @@
     """
 
-    from Numeric import array
+#    from numpy.oldnumeric import array
     import math
 
@@ -592 +591 @@
     """
 
-    from Numeric import array
     import math
 
@@ -686 +684 @@
     """
 
-    from Numeric import array
     import math
-
     from anuga.config import epsilon
-
 
     deltax = lenx/float(m)

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -11 +11 @@
 from anuga.caching import cache
 from math import pi, sqrt
-from Numeric import array, allclose
-        
+import numpy        
 
 class Mesh(General_mesh):
@@ -79 +78 @@
         """
 
-
-
-        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
-
         General_mesh.__init__(self, coordinates, triangles,
                               number_of_full_nodes=\
@@ -98 +93 @@
 
         #Allocate space for geometric quantities
-        self.centroid_coordinates = zeros((N, 2), Float)
-
-        self.radii = zeros(N, Float)
-
-        self.neighbours = zeros((N, 3), Int)
-        self.neighbour_edges = zeros((N, 3), Int)
-        self.number_of_boundaries = zeros(N, Int)
-        self.surrogate_neighbours = zeros((N, 3), Int)
+        self.centroid_coordinates = numpy.zeros((N, 2), numpy.float)
+
+        self.radii = numpy.zeros(N, numpy.float)
+
+        self.neighbours = numpy.zeros((N, 3), numpy.int)
+        self.neighbour_edges = numpy.zeros((N, 3), numpy.int)
+        self.number_of_boundaries = numpy.zeros(N, numpy.int)
+        self.surrogate_neighbours = numpy.zeros((N, 3), numpy.int)
 
         #Get x,y coordinates for all triangles and store
@@ -124 +119 @@
 
             #Compute centroid
-            centroid = array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])
+            centroid = numpy.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])
             self.centroid_coordinates[i] = centroid
 
@@ -133 +128 @@
 
                 #Midpoints
-                m0 = array([(x1 + x2)/2, (y1 + y2)/2])
-                m1 = array([(x0 + x2)/2, (y0 + y2)/2])
-                m2 = array([(x1 + x0)/2, (y1 + y0)/2])
+                m0 = numpy.array([(x1 + x2)/2, (y1 + y2)/2])
+                m1 = numpy.array([(x0 + x2)/2, (y0 + y2)/2])
+                m2 = numpy.array([(x1 + x0)/2, (y1 + y0)/2])
 
                 #The radius is the distance from the centroid of
                 #a triangle to the midpoint of the side of the triangle
                 #closest to the centroid
-                d0 = sqrt(sum( (centroid-m0)**2 ))
-                d1 = sqrt(sum( (centroid-m1)**2 ))
-                d2 = sqrt(sum( (centroid-m2)**2 ))
+                d0 = numpy.sqrt(numpy.sum( (centroid-m0)**2 ))
+                d1 = numpy.sqrt(numpy.sum( (centroid-m1)**2 ))
+                d2 = numpy.sqrt(numpy.sum( (centroid-m2)**2 ))
 
                 self.radii[i] = min(d0, d1, d2)
@@ -150 +145 @@
                 #of inscribed circle is computed
 
-                a = sqrt((x0-x1)**2+(y0-y1)**2)
-                b = sqrt((x1-x2)**2+(y1-y2)**2)
-                c = sqrt((x2-x0)**2+(y2-y0)**2)
+                a = numpy.sqrt((x0-x1)**2+(y0-y1)**2)
+                b = numpy.sqrt((x1-x2)**2+(y1-y2)**2)
+                c = numpy.sqrt((x2-x0)**2+(y2-y0)**2)
 
                 self.radii[i]=2.0*self.areas[i]/(a+b+c)
@@ -388 +383 @@
             self.element_tag is defined
         """
-        from Numeric import array, Int
 
         if tagged_elements is None:
@@ -395 +389 @@
             #Check that all keys in given boundary exist
             for tag in tagged_elements.keys():
-                tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
+                tagged_elements[tag] = numpy.array(tagged_elements[tag]).astype(numpy.int)
 
                 msg = 'Not all elements exist. '
@@ -463 +457 @@
         """
         
-        from Numeric import allclose, sqrt, array, minimum, maximum
         from anuga.utilities.numerical_tools import angle, ensure_numeric     
 
@@ -477 +470 @@
         inverse_segments = {}
         p0 = None
-        mindist = sqrt(sum((pmax-pmin)**2)) # Start value across entire mesh
+        mindist = numpy.sqrt(numpy.sum((pmax-pmin)**2)) # Start value across entire mesh
         for i, edge_id in self.boundary.keys():
             # Find vertex ids for boundary segment
@@ -487 +480 @@
             B = self.get_vertex_coordinate(i, b, absolute=True) # End
 
-
             # Take the point closest to pmin as starting point
             # Note: Could be arbitrary, but nice to have
             # a unique way of selecting
-            dist_A = sqrt(sum((A-pmin)**2))
-            dist_B = sqrt(sum((B-pmin)**2))
+            dist_A = numpy.sqrt(numpy.sum((A-pmin)**2))
+            dist_B = numpy.sqrt(numpy.sum((B-pmin)**2))
 
             # Find lower leftmost point
@@ -593 +585 @@
                 # We have reached a point already visited. 
                 
-                if allclose(p1, polygon[0]):
+                if numpy.allclose(p1, polygon[0]):
                     # If it is the initial point, the polygon is complete. 
                     
@@ -629 +621 @@
         from anuga.utilities.numerical_tools import anglediff
 
-        from Numeric import sort, allclose
-
         N = len(self)
 
@@ -672 +662 @@
 
                 # Normalise
-                l_u = sqrt(u[0]*u[0] + u[1]*u[1])
-                l_v = sqrt(v[0]*v[0] + v[1]*v[1])
+                l_u = numpy.sqrt(u[0]*u[0] + u[1]*u[1])
+                l_v = numpy.sqrt(v[0]*v[0] + v[1]*v[1])
 
                 msg = 'Normal vector in triangle %d does not have unit length' %i
-                assert allclose(l_v, 1), msg
+                assert numpy.allclose(l_v, 1), msg
 
                 x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product
@@ -730 +720 @@
 
         V = self.vertex_value_indices[:] #Take a copy
-        V = sort(V)
-        assert allclose(V, range(3*N))
+        V = numpy.sort(V)
+        assert numpy.allclose(V, range(3*N))
 
         assert sum(self.number_of_triangles_per_node) ==\
@@ -742 +732 @@
                 count[i] += 1
 
-        assert allclose(count, self.number_of_triangles_per_node)
+        assert numpy.allclose(count, self.number_of_triangles_per_node)
 
 
@@ -805 +795 @@
         """
 
-        from Numeric import arange
         from anuga.utilities.numerical_tools import histogram, create_bins
 
@@ -1141 +1130 @@
 
             # Distances from line origin to the two intersections
-            z0 = array([x0 - xi0, y0 - eta0])
-            z1 = array([x1 - xi0, y1 - eta0])              
-            d0 = sqrt(sum(z0**2)) 
-            d1 = sqrt(sum(z1**2))
+            z0 = numpy.array([x0 - xi0, y0 - eta0])
+            z1 = numpy.array([x1 - xi0, y1 - eta0])              
+            d0 = numpy.sqrt(numpy.sum(z0**2)) 
+            d1 = numpy.sqrt(numpy.sum(z1**2))
                
             if d1 < d0:
@@ -1157 +1146 @@
             # Normal direction:
             # Right hand side relative to line direction
-            vector = array([x1 - x0, y1 - y0]) # Segment vector
-            length = sqrt(sum(vector**2))      # Segment length
-            normal = array([vector[1], -vector[0]])/length
+            vector = numpy.array([x1 - x0, y1 - y0]) # Segment vector
+            length = numpy.sqrt(numpy.sum(vector**2))      # Segment length
+            normal = numpy.array([vector[1], -vector[0]])/length
 
 
@@ -1239 +1228 @@
         assert isinstance(segment, Triangle_intersection), msg
         
-        midpoint = sum(array(segment.segment))/2
+        midpoint = numpy.sum(numpy.array(segment.segment))/2
         midpoints.append(midpoint)
 

anuga_core/source/anuga/abstract_2d_finite_volumes/pmesh2domain.py

@@ -9 +9 @@
 import sys
 
+import numpy
 
 
@@ -161 +162 @@
     """
 
-    from Numeric import transpose
     from load_mesh.loadASCII import import_mesh_file
 
@@ -172 +172 @@
     volumes = mesh_dict['triangles']
     vertex_quantity_dict = {}
-    point_atts = transpose(mesh_dict['vertex_attributes'])
+    point_atts = numpy.transpose(mesh_dict['vertex_attributes'])
     point_titles  = mesh_dict['vertex_attribute_titles']
     geo_reference  = mesh_dict['geo_reference']
     if point_atts != None:
-        for quantity, value_vector in map (None, point_titles, point_atts):
+        print 'type(point_atts)=%s' % type(point_atts)
+        for quantity, value_vector in map(None, point_titles, point_atts):
             vertex_quantity_dict[quantity] = value_vector
     tag_dict = pmesh_dict_to_tag_dict(mesh_dict)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -15 +15 @@
 """
 
-from Numeric import array, zeros, Float, less, concatenate, NewAxis,\
-     argmax, argmin, allclose, take, reshape, alltrue
+import numpy
 
 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
@@ -38 +37 @@
         if vertex_values is None:
             N = len(domain) # number_of_elements
-            self.vertex_values = zeros((N, 3), Float)
+            self.vertex_values = numpy.zeros((N, 3), numpy.float)
         else:
-            self.vertex_values = array(vertex_values).astype(Float)
+            self.vertex_values = numpy.array(vertex_values).astype(numpy.float)
 
             N, V = self.vertex_values.shape
@@ -57 +56 @@
 
         # Allocate space for other quantities
-        self.centroid_values = zeros(N, Float)
-        self.edge_values = zeros((N, 3), Float)
+        self.centroid_values = numpy.zeros(N, numpy.float)
+        self.edge_values = numpy.zeros((N, 3), numpy.float)
 
         # Allocate space for Gradient
-        self.x_gradient = zeros(N, Float)
-        self.y_gradient = zeros(N, Float)
+        self.x_gradient = numpy.zeros(N, numpy.float)
+        self.y_gradient = numpy.zeros(N, numpy.float)
 
         # Allocate space for Limiter Phi
-        self.phi = zeros(N, Float)        
+        self.phi = numpy.zeros(N, numpy.float)        
 
         # Intialise centroid and edge_values
@@ -72 +71 @@
         # Allocate space for boundary values
         L = len(domain.boundary)
-        self.boundary_values = zeros(L, Float)
+        self.boundary_values = numpy.zeros(L, numpy.float)
 
         # Allocate space for updates of conserved quantities by
@@ -78 +77 @@
 
         # Allocate space for update fields
-        self.explicit_update = zeros(N, Float )
-        self.semi_implicit_update = zeros(N, Float )
-        self.centroid_backup_values = zeros(N, Float)
+        self.explicit_update = numpy.zeros(N, numpy.float )
+        self.semi_implicit_update = numpy.zeros(N, numpy.float )
+        self.centroid_backup_values = numpy.zeros(N, numpy.float)
 
         self.set_beta(1.0)
@@ -382 +381 @@
         from anuga.geospatial_data.geospatial_data import Geospatial_data
         from types import FloatType, IntType, LongType, ListType, NoneType
-        from Numeric import ArrayType
 
         # Treat special case: Polygon situation
@@ -448 +446 @@
 
         msg = 'Indices must be a list or None'
-        assert type(indices) in [ListType, NoneType, ArrayType], msg
+        assert type(indices) in [ListType, NoneType, numpy.ndarray], msg
 
 
@@ -458 +456 @@
                 self.set_values_from_constant(numeric,
                                               location, indices, verbose)
-            elif type(numeric) in [ArrayType, ListType]:
+            elif type(numeric) in [numpy.ndarray, ListType]:
                 self.set_values_from_array(numeric,
                                            location, indices, verbose)
@@ -474 +472 @@
                                                      use_cache=use_cache)
             else:
-                msg = 'Illegal type for argument numeric: %s' %str(numeric)
-                raise msg
+                msg = 'Illegal type for argument numeric: %s' % str(numeric)
+                raise TypeError, msg
 
         elif quantity is not None:
@@ -610 +608 @@
         """
 
-        from Numeric import array, Float, Int, allclose
-
-        values = array(values).astype(Float)
+        values = numpy.array(values).astype(numpy.float)
 
         if indices is not None:
-            indices = array(indices).astype(Int)
+            indices = numpy.array(indices).astype(numpy.int)
             msg = 'Number of values must match number of indices:'
             msg += ' You specified %d values and %d indices'\
@@ -640 +636 @@
 
         elif location == 'unique vertices':
-            assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
+            assert len(values.shape) == 1 or numpy.allclose(values.shape[1:], 1),\
                    'Values array must be 1d'
 
-            self.set_vertex_values(values.flat, indices=indices)
+            self.set_vertex_values(values.ravel(), indices=indices)
             
         else:
@@ -676 +672 @@
         A = q.vertex_values
 
-        from Numeric import allclose
         msg = 'Quantities are defined on different meshes. '+\
               'This might be a case for implementing interpolation '+\
               'between different meshes.'
-        assert allclose(A.shape, self.vertex_values.shape), msg
+        assert numpy.allclose(A.shape, self.vertex_values.shape), msg
 
         self.set_values(A, location='vertices',
@@ -716 +711 @@
                 indices = range(len(self))
                 
-            V = take(self.domain.get_centroid_coordinates(), indices)
+            V = numpy.take(self.domain.get_centroid_coordinates(), indices, axis=0)
+            print 'V=%s' % str(V)
             self.set_values(f(V[:,0], V[:,1]),
                             location=location,
@@ -780 +776 @@
 
 
-        points = ensure_numeric(points, Float)
-        values = ensure_numeric(values, Float)
+        points = ensure_numeric(points, numpy.float)
+        values = ensure_numeric(values, numpy.float)
 
         if location != 'vertices':
@@ -913 +909 @@
         # Always return absolute indices
         if mode is None or mode == 'max':
-            i = argmax(V)
+            i = numpy.argmax(V)
         elif mode == 'min':    
-            i = argmin(V)
+            i = numpy.argmin(V)
 
             
@@ -1118 +1114 @@
         """
         
-        from Numeric import take
-
         # FIXME (Ole): I reckon we should have the option of passing a
         #              polygon into get_values. The question becomes how
@@ -1143 +1137 @@
             raise msg
 
-        import types, Numeric
-        assert type(indices) in [types.ListType, types.NoneType,
-                                 Numeric.ArrayType],\
+        import types
+        
+        assert type(indices) in [types.ListType, types.NoneType, numpy.ndarray], \
                                  'Indices must be a list or None'
 
@@ -1151 +1145 @@
             if (indices ==  None):
                 indices = range(len(self))
-            return take(self.centroid_values,indices)
+            return numpy.take(self.centroid_values, indices, axis=0)
         elif location == 'edges':
             if (indices ==  None):
                 indices = range(len(self))
-            return take(self.edge_values,indices)
+            return numpy.take(self.edge_values, indices, axis=0)
         elif location == 'unique vertices':
             if (indices ==  None):
@@ -1178 +1172 @@
                     sum += self.vertex_values[triangle_id, vertex_id]
                 vert_values.append(sum/len(triangles))
-            return Numeric.array(vert_values)
+            return numpy.array(vert_values)
         else:
             if (indices is None):
                 indices = range(len(self))
-            return take(self.vertex_values, indices)
+            return numpy.take(self.vertex_values, indices, axis=0)
 
 
@@ -1196 +1190 @@
         """
 
-        from Numeric import array, Float
-
         # Assert that A can be converted to a Numeric array of appropriate dim
-        A = ensure_numeric(A, Float)
+        A = ensure_numeric(A, numpy.float)
 
         # print 'SHAPE A', A.shape
@@ -1273 +1265 @@
         """
 
-        from Numeric import concatenate, zeros, Float, Int, array, reshape
-
-
         if smooth is None:
             # Take default from domain
@@ -1284 +1273 @@
 
         if precision is None:
-            precision = Float
+            precision = numpy.float
             
 
@@ -1293 +1282 @@
             V = self.domain.get_triangles()
             N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
-            A = zeros(N, Float)
+            A = numpy.zeros(N, numpy.float)
             points = self.domain.get_nodes()            
             
@@ -1343 +1332 @@
             V = self.domain.get_disconnected_triangles()
             points = self.domain.get_vertex_coordinates()
-            A = self.vertex_values.flat.astype(precision)
+            A = self.vertex_values.ravel().astype(precision)
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -11 +11 @@
 
 #include "Python.h"
-#include "Numeric/arrayobject.h"
+#include "numpy/arrayobject.h"
 #include "math.h"
 

anuga_core/source/anuga/abstract_2d_finite_volumes/region.py

@@ -8 +8 @@
 # FIXME (DSG-DSG) add better comments
 
-from Numeric import average
+import numpy
+
+
 class Region:
     """Base class for modifying quantities based on a region.
@@ -101 +103 @@
                 values = Q.get_values(indices=self.build_indices(elements, domain),
                                       location=self.location)
-                av = average(values)
+                av = numpy.average(values)
                 if self.location == "vertices":
-                    av = average(av)
+                    av = numpy.average(av)
                 new_values = av + self.X    
             else:

anuga_core/source/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py

@@ -18 +18 @@
 
 from mesh_factory import rectangular
-from Numeric import array
-    
+
 
 ######################
@@ -78 +77 @@
 domain.set_quantity('stage', Z)
 
-from Numeric import allclose
-
 #Evolve
 for t in domain.evolve(yieldstep = 0.1, finaltime = 30):

anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

@@ -6 +6 @@
 from domain import *
 from anuga.config import epsilon
-from Numeric import allclose, array, ones, Float
+import numpy
 
 
@@ -64 +64 @@
 
 
-        assert domain.get_conserved_quantities(0, edge=1) == 0.
+        assert numpy.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
 
 
@@ -95 +95 @@
         #Centroids
         q = domain.get_conserved_quantities(0)
-        assert allclose(q, [2., 2., 0.])
+        assert numpy.allclose(q, [2., 2., 0.])
 
         q = domain.get_conserved_quantities(1)
-        assert allclose(q, [5., 5., 0.])
+        assert numpy.allclose(q, [5., 5., 0.])
 
         q = domain.get_conserved_quantities(2)
-        assert allclose(q, [3., 3., 0.])
+        assert numpy.allclose(q, [3., 3., 0.])
 
         q = domain.get_conserved_quantities(3)
-        assert allclose(q, [0., 0., 0.])
+        assert numpy.allclose(q, [0., 0., 0.])
 
 
         #Edges
         q = domain.get_conserved_quantities(0, edge=0)
-        assert allclose(q, [2.5, 2.5, 0.])
+        assert numpy.allclose(q, [2.5, 2.5, 0.])
         q = domain.get_conserved_quantities(0, edge=1)
-        assert allclose(q, [2., 2., 0.])
+        assert numpy.allclose(q, [2., 2., 0.])
         q = domain.get_conserved_quantities(0, edge=2)
-        assert allclose(q, [1.5, 1.5, 0.])
+        assert numpy.allclose(q, [1.5, 1.5, 0.])
 
         for i in range(3):
             q = domain.get_conserved_quantities(1, edge=i)
-            assert allclose(q, [5, 5, 0.])
+            assert numpy.allclose(q, [5, 5, 0.])
 
 
         q = domain.get_conserved_quantities(2, edge=0)
-        assert allclose(q, [4.5, 4.5, 0.])
+        assert numpy.allclose(q, [4.5, 4.5, 0.])
         q = domain.get_conserved_quantities(2, edge=1)
-        assert allclose(q, [4.5, 4.5, 0.])
+        assert numpy.allclose(q, [4.5, 4.5, 0.])
         q = domain.get_conserved_quantities(2, edge=2)
-        assert allclose(q, [0., 0., 0.])
+        assert numpy.allclose(q, [0., 0., 0.])
 
 
         q = domain.get_conserved_quantities(3, edge=0)
-        assert allclose(q, [3., 3., 0.])
+        assert numpy.allclose(q, [3., 3., 0.])
         q = domain.get_conserved_quantities(3, edge=1)
-        assert allclose(q, [-1.5, -1.5, 0.])
+        assert numpy.allclose(q, [-1.5, -1.5, 0.])
         q = domain.get_conserved_quantities(3, edge=2)
-        assert allclose(q, [-1.5, -1.5, 0.])
+        assert numpy.allclose(q, [-1.5, -1.5, 0.])
 
 
@@ -179 +179 @@
         Q = domain.create_quantity_from_expression(expression)
 
-        assert allclose(Q.vertex_values, [[2,3,4], [6,6,6],
-                                      [1,1,10], [-5, 4, 4]])
+        assert numpy.allclose(Q.vertex_values, [[2,3,4], [6,6,6],
+                                                [1,1,10], [-5, 4, 4]])
 
         expression = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
@@ -188 +188 @@
         Y = domain.quantities['ymomentum'].vertex_values
 
-        assert allclose(Q.vertex_values, (X**2 + Y**2)**0.5)
+        assert numpy.allclose(Q.vertex_values, (X**2 + Y**2)**0.5)
 
 
@@ -314 +314 @@
         Q = domain.quantities['depth']
 
-        assert allclose(Q.vertex_values, [[2,3,4], [6,6,6],
-                                      [1,1,10], [-5, 4, 4]])
+        assert numpy.allclose(Q.vertex_values, [[2,3,4], [6,6,6],
+                                                [1,1,10], [-5, 4, 4]])
 
 
@@ -382 +382 @@
         domain.check_integrity()
 
-        assert allclose(domain.neighbours, [[-1,-2,-3]])
+        assert numpy.allclose(domain.neighbours, [[-1,-2,-3]])
 
 
@@ -471 +471 @@
                                       [0,0,9], [-6, 3, 3]])
 
-        assert allclose( domain.quantities['stage'].centroid_values,
-                         [2,5,3,0] )
+        assert numpy.allclose( domain.quantities['stage'].centroid_values,
+                               [2,5,3,0] )
 
         domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
@@ -484 +484 @@
 
         #First order extrapolation
-        assert allclose( domain.quantities['stage'].vertex_values,
-                         [[ 2.,  2.,  2.],
-                          [ 5.,  5.,  5.],
-                          [ 3.,  3.,  3.],
-                          [ 0.,  0.,  0.]])
+        assert numpy.allclose( domain.quantities['stage'].vertex_values,
+                               [[ 2.,  2.,  2.],
+                                [ 5.,  5.,  5.],
+                                [ 3.,  3.,  3.],
+                                [ 0.,  0.,  0.]])
 
 
@@ -527 +527 @@
 
         for name in domain.conserved_quantities:
-            domain.quantities[name].explicit_update = array([4.,3.,2.,1.])
-            domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.])
+            domain.quantities[name].explicit_update = numpy.array([4.,3.,2.,1.])
+            domain.quantities[name].semi_implicit_update = numpy.array([1.,1.,1.,1.])
 
 
@@ -535 +535 @@
         domain.update_conserved_quantities()
 
-        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
-        denom = ones(4, Float)-domain.timestep*sem
+        sem = numpy.array([1.,1.,1.,1.])/numpy.array([1, 2, 3, 4])
+        denom = numpy.ones(4, numpy.float)-domain.timestep*sem
 
 #        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
 #        x /= denom
 
-        x = array([1., 2., 3., 4.])
+        x = numpy.array([1., 2., 3., 4.])
         x /= denom
-        x += domain.timestep*array( [4,3,2,1] )
+        x += domain.timestep*numpy.array( [4,3,2,1] )
 
         for name in domain.conserved_quantities:
-            assert allclose(domain.quantities[name].centroid_values, x)
+            assert numpy.allclose(domain.quantities[name].centroid_values, x)
 
 
@@ -578 +578 @@
                                       [0,0,9], [-6, 3, 3]])
 
-        assert allclose( domain.quantities['stage'].centroid_values,
-                         [2,5,3,0] )
+        assert numpy.allclose( domain.quantities['stage'].centroid_values,
+                               [2,5,3,0] )
 
         domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
@@ -591 +591 @@
 
         #First order extrapolation
-        assert allclose( domain.quantities['stage'].vertex_values,
-                         [[ 2.,  2.,  2.],
-                          [ 5.,  5.,  5.],
-                          [ 3.,  3.,  3.],
-                          [ 0.,  0.,  0.]])
+        assert numpy.allclose( domain.quantities['stage'].vertex_values,
+                               [[ 2.,  2.,  2.],
+                                [ 5.,  5.,  5.],
+                                [ 3.,  3.,  3.],
+                                [ 0.,  0.,  0.]])
 
         domain.build_tagged_elements_dictionary({'mound':[0,1]})
@@ -610 +610 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -627 +626 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -645 +643 @@
         domain.set_region([set_bottom_friction, set_top_friction])
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),\
-                        [[ 0.09,  0.09,  0.09],
-                         [ 0.09,  0.09,  0.09],
-                         [ 0.07,  0.07,  0.07],
-                         [ 0.07,  0.07,  0.07],
-                         [ 1.0,  1.0,  1.0],
-                         [ 1.0,  1.0,  1.0]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),\
+                              [[ 0.09,  0.09,  0.09],
+                               [ 0.09,  0.09,  0.09],
+                               [ 0.07,  0.07,  0.07],
+                               [ 0.07,  0.07,  0.07],
+                               [ 1.0,  1.0,  1.0],
+                               [ 1.0,  1.0,  1.0]])
 
         domain.set_region([set_all_friction])
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),
-                        [[ 10.09, 10.09, 10.09],
-                         [ 10.09, 10.09, 10.09],
-                         [ 10.07, 10.07, 10.07],
-                         [ 10.07, 10.07, 10.07],
-                         [ 11.0,  11.0,  11.0],
-                         [ 11.0,  11.0,  11.0]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),
+                              [[ 10.09, 10.09, 10.09],
+                               [ 10.09, 10.09, 10.09],
+                               [ 10.07, 10.07, 10.07],
+                               [ 10.07, 10.07, 10.07],
+                               [ 11.0,  11.0,  11.0],
+                               [ 11.0,  11.0,  11.0]])
 
 
@@ -670 +668 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -690 +687 @@
         
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),\
-                        [[ 0.09,  0.09,  0.09],
-                         [ 0.09,  0.09,  0.09],
-                         [ 0.07,  0.07,  0.07],
-                         [ 0.07,  0.07,  0.07],
-                         [ 1.0,  1.0,  1.0],
-                         [ 1.0,  1.0,  1.0]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),
+                              [[ 0.09,  0.09,  0.09],
+                               [ 0.09,  0.09,  0.09],
+                               [ 0.07,  0.07,  0.07],
+                               [ 0.07,  0.07,  0.07],
+                               [ 1.0,  1.0,  1.0],
+                               [ 1.0,  1.0,  1.0]])
         
         domain.set_region([set_bottom_friction, set_top_friction])
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),\
-                        [[ 0.09,  0.09,  0.09],
-                         [ 0.09,  0.09,  0.09],
-                         [ 0.07,  0.07,  0.07],
-                         [ 0.07,  0.07,  0.07],
-                         [ 1.0,  1.0,  1.0],
-                         [ 1.0,  1.0,  1.0]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),
+                              [[ 0.09,  0.09,  0.09],
+                               [ 0.09,  0.09,  0.09],
+                               [ 0.07,  0.07,  0.07],
+                               [ 0.07,  0.07,  0.07],
+                               [ 1.0,  1.0,  1.0],
+                               [ 1.0,  1.0,  1.0]])
 
         domain.set_region([set_all_friction])
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),
-                        [[ 10.09, 10.09, 10.09],
-                         [ 10.09, 10.09, 10.09],
-                         [ 10.07, 10.07, 10.07],
-                         [ 10.07, 10.07, 10.07],
-                         [ 11.0,  11.0,  11.0],
-                         [ 11.0,  11.0,  11.0]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),
+                              [[ 10.09, 10.09, 10.09],
+                               [ 10.09, 10.09, 10.09],
+                               [ 10.07, 10.07, 10.07],
+                               [ 10.07, 10.07, 10.07],
+                               [ 11.0,  11.0,  11.0],
+                               [ 11.0,  11.0,  11.0]])
 
 #-------------------------------------------------------------

anuga_core/source/anuga/abstract_2d_finite_volumes/test_ermapper.py

@@ -5 +5 @@
 
 import ermapper_grids
-import Numeric
+import numpy
 from os import remove
 
@@ -19 +19 @@
         data_filename = 'test_write_ermapper_grid'
         
-        original_grid = Numeric.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]])
+        original_grid = numpy.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]])
 
         # Check that the function works when passing the filename without
@@ -25 +25 @@
         ermapper_grids.write_ermapper_grid(data_filename, original_grid)
         new_grid = ermapper_grids.read_ermapper_grid(data_filename)
-        assert Numeric.allclose(original_grid, new_grid)
+        assert numpy.allclose(original_grid, new_grid)
 
         # Check that the function works when passing the filename with
@@ -31 +31 @@
         ermapper_grids.write_ermapper_grid(header_filename, original_grid)
         new_grid = ermapper_grids.read_ermapper_grid(header_filename)
-        assert Numeric.allclose(original_grid, new_grid)
+        assert numpy.allclose(original_grid, new_grid)
 
         # Clean up created files
@@ -40 +40 @@
         # Setup test data
         filename = 'test_write_ermapper_grid'
-        original_grid = Numeric.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])
+        original_grid = numpy.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])
 
         # Write test data
-        ermapper_grids.write_ermapper_data(original_grid, filename, Numeric.Float64)
+        ermapper_grids.write_ermapper_data(original_grid, filename, numpy.float64)
 
         # Read in the test data
-        new_grid = ermapper_grids.read_ermapper_data(filename, Numeric.Float64)
+        new_grid = ermapper_grids.read_ermapper_data(filename, numpy.float64)
 
         # Check that the test data that has been read in matches the original data
-        assert Numeric.allclose(original_grid, new_grid)
+        assert numpy.allclose(original_grid, new_grid)
 
         # Clean up created files
@@ -57 +57 @@
         # Setup test data
         filename = 'test_write_ermapper_grid'
-        original_grid = Numeric.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])
+        original_grid = numpy.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])
 
         # Write test data
@@ -66 +66 @@
 
         # Check that the test data that has been read in matches the original data
-        assert Numeric.allclose(original_grid, new_grid)
+        assert numpy.allclose(original_grid, new_grid)
 
         # Clean up created files
@@ -75 +75 @@
 
         # setup test data
-        original_grid = Numeric.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]])
+        original_grid = numpy.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]])
         # Write test data
         ermapper_grids.write_ermapper_data(original_grid, data_filename)

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -7 +7 @@
 
 from anuga.config import epsilon
-from Numeric import allclose, array, ones, Float
+import numpy
 from general_mesh import General_mesh
 from anuga.coordinate_transforms.geo_reference import Geo_reference
@@ -23 +23 @@
     def test_get_vertex_coordinates(self):
         from mesh_factory import rectangular
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -30 +29 @@
 
 
-        assert allclose(domain.get_nodes(), nodes)
+        assert numpy.allclose(domain.get_nodes(), nodes)
 
 
@@ -41 +40 @@
             for j in range(3):
                 k = triangles[i,j]  #Index of vertex j in triangle i
-                assert allclose(V[3*i+j,:], nodes[k])
+                assert numpy.allclose(V[3*i+j,:], nodes[k])
 
     def test_get_vertex_coordinates_with_geo_ref(self):
@@ -55 +54 @@
         f = [4.0, 0.0]
 
-        nodes = array([a, b, c, d, e, f])
+        nodes = numpy.array([a, b, c, d, e, f])
 
         nodes_absolute = geo.get_absolute(nodes)
         
         #bac, bce, ecf, dbe, daf, dae
-        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        triangles = numpy.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
         domain = General_mesh(nodes, triangles,
                        geo_reference = geo)
         verts = domain.get_vertex_coordinates(triangle_id=0)        
-        self.assert_(allclose(array([b,a,c]), verts)) 
+        self.assert_(numpy.allclose(numpy.array([b,a,c]), verts)) 
         verts = domain.get_vertex_coordinates(triangle_id=0)       
-        self.assert_(allclose(array([b,a,c]), verts))
+        self.assert_(numpy.allclose(numpy.array([b,a,c]), verts))
         verts = domain.get_vertex_coordinates(triangle_id=0,
                                               absolute=True)       
-        self.assert_(allclose(array([nodes_absolute[1],
-                                     nodes_absolute[0],
-                                     nodes_absolute[2]]), verts))
+        self.assert_(numpy.allclose(numpy.array([nodes_absolute[1],
+                                                 nodes_absolute[0],
+                                                 nodes_absolute[2]]), verts))
         verts = domain.get_vertex_coordinates(triangle_id=0,
                                               absolute=True)       
-        self.assert_(allclose(array([nodes_absolute[1],
-                                     nodes_absolute[0],
-                                     nodes_absolute[2]]), verts))
+        self.assert_(numpy.allclose(numpy.array([nodes_absolute[1],
+                                                 nodes_absolute[0],
+                                                 nodes_absolute[2]]), verts))
         
         
@@ -86 +85 @@
         """
         from mesh_factory import rectangular
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -93 +91 @@
 
 
-        assert allclose(domain.get_nodes(), nodes)
+        assert numpy.allclose(domain.get_nodes(), nodes)
 
 
@@ -104 +102 @@
             for j in range(3):
                 k = triangles[i,j]  #Index of vertex j in triangle i
-                assert allclose(V[j,:], nodes[k])
+                assert numpy.allclose(V[j,:], nodes[k])
 
 
@@ -114 +112 @@
         """
         from mesh_factory import rectangular
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -121 +118 @@
 
         value = [7]
-        assert allclose(domain.get_triangles(), triangles)
-        assert allclose(domain.get_triangles([0,4]),
-                        [triangles[0], triangles[4]])
+        assert numpy.allclose(domain.get_triangles(), triangles)
+        assert numpy.allclose(domain.get_triangles([0,4]),
+                              [triangles[0], triangles[4]])
         
 
@@ -130 +127 @@
         """
         from mesh_factory import rectangular
-        from Numeric import zeros, Float, array
-
-        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]
-
-        nodes = array([a, b, c, d, e, f])
+
+        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]
+
+        nodes = numpy.array([a, b, c, d, e, f])
         #bac, bce, ecf, dbe, daf, dae
-        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        triangles = numpy.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
         domain1 = General_mesh(nodes, triangles)
@@ -162 +158 @@
             #print count
             #
-            assert allclose(count, domain.number_of_triangles_per_node)
+            assert numpy.allclose(count, domain.number_of_triangles_per_node)
             
             # Check indices
@@ -196 +192 @@
         f = [4.0, 0.0]
 
-        nodes = array([a, b, c, d, e, f])
+        nodes = numpy.array([a, b, c, d, e, f])
         #bac, bce, ecf, dbe, daf, dae
-        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        triangles = numpy.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
         domain = General_mesh(nodes, triangles)
@@ -204 +200 @@
         # One node
         L = domain.get_triangles_and_vertices_per_node(node=2)
-        assert allclose(L[0], [0, 2])
-        assert allclose(L[1], [1, 1])
-        assert allclose(L[2], [2, 1])
+        print 'L=%s' % str(L)
+        assert numpy.allclose(L[0], [0, 2])
+        assert numpy.allclose(L[1], [1, 1])
+        assert numpy.allclose(L[2], [2, 1])
 
         # All nodes
@@ -213 +210 @@
         for i, Lref in enumerate(ALL):
             L = domain.get_triangles_and_vertices_per_node(node=i)
-            assert allclose(L, Lref)
+            assert numpy.allclose(L, Lref)
             
 
@@ -224 +221 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -239 +235 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -273 +268 @@
         f = [4.0, 0.0]
 
-        nodes = array([a, b, c, d, e, f])
+        nodes = numpy.array([a, b, c, d, e, f])
 
         nodes_absolute = geo.get_absolute(nodes)
         
         #bac, bce, ecf, dbe, daf, dae
-        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        triangles = numpy.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
         domain = General_mesh(nodes, triangles,
                        geo_reference = geo)
         node = domain.get_node(2)        
-        self.assertEqual(c, node)
+        self.assertTrue(numpy.alltrue(c == node))
         
         node = domain.get_node(2, absolute=True)     
-        self.assertEqual(nodes_absolute[2], node)
+        self.assertTrue(numpy.alltrue(nodes_absolute[2] == node))
         
         node = domain.get_node(2, absolute=True)     
-        self.assertEqual(nodes_absolute[2], node)
+        self.assertTrue(numpy.alltrue(nodes_absolute[2] == node))
         
 
@@ -310 +305 @@
         f = [4.0, 0.0]
 
-        nodes = array([a, b, c, d, e, f])
+        nodes = numpy.array([a, b, c, d, e, f])
 
         nodes_absolute = geo.get_absolute(nodes)
         
         # max index is 5, use 5, expect success
-        triangles = array([[1,5,2], [1,2,4], [4,2,5], [3,1,4]])
+        triangles = numpy.array([[1,5,2], [1,2,4], [4,2,5], [3,1,4]])
         General_mesh(nodes, triangles, geo_reference=geo)
         
         # max index is 5, use 6, expect assert failure
-        triangles = array([[1,6,2], [1,2,4], [4,2,5], [3,1,4]])
+        triangles = numpy.array([[1,6,2], [1,2,4], [4,2,5], [3,1,4]])
         self.failUnlessRaises(AssertionError, General_mesh,
                               nodes, triangles, geo_reference=geo)
         
         # max index is 5, use 10, expect assert failure
-        triangles = array([[1,10,2], [1,2,4], [4,2,5], [3,1,4]])
+        triangles = numpy.array([[1,10,2], [1,2,4], [4,2,5], [3,1,4]])
         self.failUnlessRaises(AssertionError, General_mesh,
                               nodes, triangles, geo_reference=geo)
@@ -333 +328 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_General_Mesh,'test') 
-    #suite = unittest.makeSuite(Test_General_Mesh,'test_get_node')    
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

@@ -6 +6 @@
 from generic_boundary_conditions import *
 from anuga.config import epsilon
-from Numeric import allclose, array
+import numpy
 
 
@@ -44 +44 @@
 
         q = Bd.evaluate()
-        assert allclose(q, x)
+        assert numpy.allclose(q, x)
 
 
@@ -98 +98 @@
         q = T.evaluate(0, 2)  #Vol=0, edge=2
 
-        assert allclose(q, [1.5, 2.5])
+        assert numpy.allclose(q, [1.5, 2.5])
 
 
@@ -175 +175 @@
 
         #Check that midpoint coordinates at boundary are correctly computed
-        assert allclose( F.midpoint_coordinates,
-                         [[1.0, 0.0], [0.0, 1.0], [3.0, 0.0],
-                          [3.0, 1.0], [1.0, 3.0], [0.0, 3.0]])
-
-        #assert allclose(F.midpoint_coordinates[(3,2)], [0.0, 3.0])
-        #assert allclose(F.midpoint_coordinates[(3,1)], [1.0, 3.0])
-        #assert allclose(F.midpoint_coordinates[(0,2)], [0.0, 1.0])
-        #assert allclose(F.midpoint_coordinates[(0,0)], [1.0, 0.0])
-        #assert allclose(F.midpoint_coordinates[(2,0)], [3.0, 0.0])
-        #assert allclose(F.midpoint_coordinates[(2,1)], [3.0, 1.0])
+        assert numpy.allclose( F.midpoint_coordinates,
+                               [[1.0, 0.0], [0.0, 1.0], [3.0, 0.0],
+                                [3.0, 1.0], [1.0, 3.0], [0.0, 3.0]])
+
+        #assert numpy.allclose(F.midpoint_coordinates[(3,2)], [0.0, 3.0])
+        #assert numpy.allclose(F.midpoint_coordinates[(3,1)], [1.0, 3.0])
+        #assert numpy.allclose(F.midpoint_coordinates[(0,2)], [0.0, 1.0])
+        #assert numpy.allclose(F.midpoint_coordinates[(0,0)], [1.0, 0.0])
+        #assert numpy.allclose(F.midpoint_coordinates[(2,0)], [3.0, 0.0])
+        #assert numpy.allclose(F.midpoint_coordinates[(2,1)], [3.0, 1.0])
 
 
@@ -193 +193 @@
         domain.time = 5*30/2  #A quarter way through first step
         q = F.evaluate()
-        assert allclose(q, [1.0/4, sin(2*pi/10)/4])
+        assert numpy.allclose(q, [1.0/4, sin(2*pi/10)/4])
 
 
         domain.time = 2.5*5*60  #Half way between steps 2 and 3
         q = F.evaluate()
-        assert allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2])
+        assert numpy.allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2])
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/test_ghost.py

@@ -6 +6 @@
 from domain import *
 from anuga.config import epsilon
-from Numeric import allclose, array, ones, Float
-
-
+import numpy
 
 
@@ -43 +41 @@
 
 
-        assert domain.get_conserved_quantities(0, edge=1) == 0.
+        assert numpy.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -1 +1 @@
 #!/usr/bin/env python
-
-
 
 #FIXME: Seperate the tests for mesh and general_mesh
@@ -13 +11 @@
 from mesh_factory import rectangular
 from anuga.config import epsilon
-from Numeric import allclose, array, Int
+import numpy
 
 from anuga.coordinate_transforms.geo_reference import Geo_reference
@@ -20 +18 @@
 
 def distance(x, y):
-    return sqrt( sum( (array(x)-array(y))**2 ))
+    return sqrt( sum( (numpy.array(x)-numpy.array(y))**2 ))
 
 class Test_Mesh(unittest.TestCase):
@@ -63 +61 @@
         #Normals
         normals = mesh.get_normals()
-        assert allclose(normals[0, 0:2], [3.0/5, 4.0/5])
-        assert allclose(normals[0, 2:4], [-1.0, 0.0])
-        assert allclose(normals[0, 4:6], [0.0, -1.0])
-
-        assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])
-        assert allclose(mesh.get_normal(0,1), [-1.0, 0.0])
-        assert allclose(mesh.get_normal(0,2), [0.0, -1.0])
+        assert numpy.allclose(normals[0, 0:2], [3.0/5, 4.0/5])
+        assert numpy.allclose(normals[0, 2:4], [-1.0, 0.0])
+        assert numpy.allclose(normals[0, 4:6], [0.0, -1.0])
+
+        assert numpy.allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])
+        assert numpy.allclose(mesh.get_normal(0,1), [-1.0, 0.0])
+        assert numpy.allclose(mesh.get_normal(0,2), [0.0, -1.0])
 
         #Edge lengths
-        assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])
+        assert numpy.allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])
 
 
@@ -81 +79 @@
 
         V = mesh.get_vertex_coordinates()
-        assert allclose(V, [ [0.0, 0.0],
-                             [4.0, 0.0],
-                             [0.0, 3.0] ])
+        assert numpy.allclose(V, [ [0.0, 0.0],
+                                   [4.0, 0.0],
+                                   [0.0, 3.0] ])
 
         V0 = mesh.get_vertex_coordinate(0, 0)
-        assert allclose(V0, [0.0, 0.0])
+        assert numpy.allclose(V0, [0.0, 0.0])
 
         V1 = mesh.get_vertex_coordinate(0, 1)
-        assert allclose(V1, [4.0, 0.0])
+        assert numpy.allclose(V1, [4.0, 0.0])
 
         V2 = mesh.get_vertex_coordinate(0, 2)
-        assert allclose(V2, [0.0, 3.0])
+        assert numpy.allclose(V2, [0.0, 3.0])
 
 
@@ -237 +235 @@
 
         mesh = Mesh(points, vertices,use_inscribed_circle=False)
-        assert allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work'
+        assert numpy.allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work'
 
         mesh = Mesh(points, vertices,use_inscribed_circle=True)
-        assert allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work'
+        assert numpy.allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work'
 
     def test_inscribed_circle_rightangle_triangle(self):
@@ -252 +250 @@
 
         mesh = Mesh(points, vertices,use_inscribed_circle=False)
-        assert allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work'
+        assert numpy.allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work'
 
         mesh = Mesh(points, vertices,use_inscribed_circle=True)
-        assert allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work'
+        assert numpy.allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work'
 
 
@@ -269 +267 @@
         assert mesh.areas[0] == 2.0
 
-        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
+        assert numpy.allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
 
 
@@ -303 +301 @@
         assert mesh.edgelengths[1,2] == sqrt(8.0)
 
-        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
-        assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])
-        assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])
-        assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])
+        assert numpy.allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
+        assert numpy.allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])
+        assert numpy.allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])
+        assert numpy.allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])
 
     def test_mesh_and_neighbours(self):
@@ -697 +695 @@
         get values based on triangle lists.
         """
+
         from mesh_factory import rectangular
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -716 +714 @@
     def test_boundary_polygon(self):
         from mesh_factory import rectangular
-        #from mesh import Mesh
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -727 +723 @@
 
         assert len(P) == 8
-        assert allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],
-                            [1.0, 0.5], [1.0, 1.0], [0.5, 1.0],
-                            [0.0, 1.0], [0.0, 0.5]])
+        assert numpy.allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],
+                                  [1.0, 0.5], [1.0, 1.0], [0.5, 1.0],
+                                  [0.0, 1.0], [0.0, 0.5]])
         for p in points:
             #print p, P
@@ -736 +732 @@
 
     def test_boundary_polygon_II(self):
-        from Numeric import zeros, Float
-        
-
         #Points
         a = [0.0, 0.0] #0
@@ -764 +757 @@
 
         assert len(P) == 8
-        assert allclose(P, [a, d, f, i, h, e, c, b])
+        assert numpy.allclose(P, [a, d, f, i, h, e, c, b])
 
         for p in points:
@@ -774 +767 @@
         """Same as II but vertices ordered differently
         """
-
-        from Numeric import zeros, Float
-
 
         #Points
@@ -804 +794 @@
 
         assert len(P) == 8
-        assert allclose(P, [a, d, f, i, h, e, c, b])
+        assert numpy.allclose(P, [a, d, f, i, h, e, c, b])
 
         for p in points:
@@ -815 +805 @@
         is partitioned using pymetis.
         """
-
-        from Numeric import zeros, Float
-
 
         #Points
@@ -847 +834 @@
         
         # Note that point e appears twice!
-        assert allclose(P, [a, d, f, e, g, h, e, c, b])
+        assert numpy.allclose(P, [a, d, f, e, g, h, e, c, b])
 
         for p in points:
@@ -864 +851 @@
         """
 
-        from Numeric import zeros, Float
         from mesh_factory import rectangular        
 
@@ -907 +893 @@
         
         """
-        from Numeric import zeros, Float
-        
 
         #Points
@@ -955 +939 @@
         
         assert len(P) == 8
-        assert allclose(P, [a, d, f, i, h, e, c, b])
-        assert allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)])
+        assert numpy.allclose(P, [a, d, f, i, h, e, c, b])
+        assert numpy.allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)])
         
 
@@ -998 +982 @@
                   [  52341.70703125,  38563.39453125]]
 
-        ##points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation        
+        ##points = ensure_numeric(points, int)/1000  # Simplify for ease of interpretation        
 
         triangles = [[19, 0,15],
@@ -1101 +1085 @@
                   [  35406.3359375 ,  79332.9140625 ]]
 
-        scaled_points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation
+        scaled_points = ensure_numeric(points, int)/1000  # Simplify for ease of interpretation
 
         triangles = [[ 0, 1, 2],
@@ -1142 +1126 @@
             assert is_inside_polygon(p, P)            
 
-        assert allclose(P, Pref)    
+        assert numpy.allclose(P, Pref)    
 
     def test_lone_vertices(self):
@@ -1194 +1178 @@
         boundary_polygon = mesh.get_boundary_polygon()
 
-        assert allclose(absolute_points, boundary_polygon)
+        assert numpy.allclose(absolute_points, boundary_polygon)
 
     def test_get_triangle_containing_point(self):
@@ -1241 +1225 @@
 
         neighbours = mesh.get_triangle_neighbours(0)
-        assert allclose(neighbours, [-1,1,-1])
+        assert numpy.allclose(neighbours, [-1,1,-1])
         neighbours = mesh.get_triangle_neighbours(-10)
         assert neighbours == []
@@ -1290 +1274 @@
             for x in L:
                 if x.triangle_id % 2 == 0:
-                    assert allclose(x.length, ceiling-y_line)
+                    assert numpy.allclose(x.length, ceiling-y_line)
                 else:
-                    assert allclose(x.length, y_line-floor)                
+                    assert numpy.allclose(x.length, y_line-floor)                
 
                 
-                assert allclose(x.normal, [0,-1])
-
-                assert allclose(x.segment[1][0], x.segment[0][0] + x.length)
-                assert allclose(x.segment[0][1], y_line)
-                assert allclose(x.segment[1][1], y_line)                
+                assert numpy.allclose(x.normal, [0,-1])
+
+                assert numpy.allclose(x.segment[1][0], x.segment[0][0] + x.length)
+                assert numpy.allclose(x.segment[0][1], y_line)
+                assert numpy.allclose(x.segment[1][1], y_line)                
 
                 assert x.triangle_id in intersected_triangles
@@ -1306 +1290 @@
 
             msg = 'Segments do not add up'
-            assert allclose(total_length, 2), msg
+            assert numpy.allclose(total_length, 2), msg
             
 
@@ -1341 +1325 @@
         total_length = 0
         for x in L:
-            assert allclose(x.length, 1.0)
-            assert allclose(x.normal, [0,-1])
-
-            assert allclose(x.segment[1][0], x.segment[0][0] + x.length)
-            assert allclose(x.segment[0][1], y_line)
-            assert allclose(x.segment[1][1], y_line)                            
+            assert numpy.allclose(x.length, 1.0)
+            assert numpy.allclose(x.normal, [0,-1])
+
+            assert numpy.allclose(x.segment[1][0], x.segment[0][0] + x.length)
+            assert numpy.allclose(x.segment[0][1], y_line)
+            assert numpy.allclose(x.segment[1][1], y_line)                            
 
 
@@ -1355 +1339 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 2), msg
+        assert numpy.allclose(total_length, 2), msg
         
 
@@ -1395 +1379 @@
         for x in L:
             if x.triangle_id == 1:
-                assert allclose(x.length, 1)        
-                assert allclose(x.normal, [0, -1])
+                assert numpy.allclose(x.length, 1)        
+                assert numpy.allclose(x.normal, [0, -1])
                 
             if x.triangle_id == 5:
-                assert allclose(x.length, 0.5)
-                assert allclose(x.normal, [0, -1])
+                assert numpy.allclose(x.length, 0.5)
+                assert numpy.allclose(x.normal, [0, -1])
 
 
@@ -1408 +1392 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 1.5), msg            
+        assert numpy.allclose(total_length, 1.5), msg            
 
 
@@ -1442 +1426 @@
         total_length = 0
         for i, x in enumerate(L): 
-            assert allclose(x.length, s2)
-            assert allclose(x.normal, [-s2, -s2])
-            assert allclose(sum(x.normal**2), 1)
+            assert numpy.allclose(x.length, s2)
+            assert numpy.allclose(x.normal, [-s2, -s2])
+            assert numpy.allclose(sum(x.normal**2), 1)
             
             assert x.triangle_id in intersected_triangles
@@ -1451 +1435 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 4*s2), msg
+        assert numpy.allclose(total_length, 4*s2), msg
 
 
@@ -1466 +1450 @@
         total_length = 0
         for i, x in enumerate(L): 
-            assert allclose(x.length, s2)
-            assert allclose(x.normal, [s2, s2])
-            assert allclose(sum(x.normal**2), 1)
+            assert numpy.allclose(x.length, s2)
+            assert numpy.allclose(x.normal, [s2, s2])
+            assert numpy.allclose(sum(x.normal**2), 1)
             
             assert x.triangle_id in intersected_triangles
@@ -1475 +1459 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 4*s2), msg                    
+        assert numpy.allclose(total_length, 4*s2), msg                    
 
 
@@ -1491 +1475 @@
         total_length = 0
         for i, x in enumerate(L): 
-            assert allclose(x.length, 2*s2)
-            assert allclose(x.normal, [-s2, s2])
-            assert allclose(sum(x.normal**2), 1)
+            assert numpy.allclose(x.length, 2*s2)
+            assert numpy.allclose(x.normal, [-s2, s2])
+            assert numpy.allclose(sum(x.normal**2), 1)
             
             assert x.triangle_id in intersected_triangles
@@ -1500 +1484 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 4*s2), msg                        
+        assert numpy.allclose(total_length, 4*s2), msg                        
 
 
@@ -1515 +1499 @@
         total_length = 0
         for i, x in enumerate(L): 
-            assert allclose(x.length, 2*s2)
-            assert allclose(x.normal, [s2, -s2])
-            assert allclose(sum(x.normal**2), 1)
+            assert numpy.allclose(x.length, 2*s2)
+            assert numpy.allclose(x.normal, [s2, -s2])
+            assert numpy.allclose(sum(x.normal**2), 1)
             
             assert x.triangle_id in intersected_triangles
@@ -1524 +1508 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 4*s2), msg                        
+        assert numpy.allclose(total_length, 4*s2), msg                        
 
 
@@ -1540 +1524 @@
         total_length = 0
         for i, x in enumerate(L): 
-            assert allclose(x.length, s2)
-            assert allclose(x.normal, [-s2, -s2])
-            assert allclose(sum(x.normal**2), 1)
+            assert numpy.allclose(x.length, s2)
+            assert numpy.allclose(x.normal, [-s2, -s2])
+            assert numpy.allclose(sum(x.normal**2), 1)
             
             assert x.triangle_id in intersected_triangles
@@ -1549 +1533 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 2*s2), msg
+        assert numpy.allclose(total_length, 2*s2), msg
 
 
@@ -1562 +1546 @@
         total_length = 0
         for i, x in enumerate(L): 
-            assert allclose(x.normal, [-s2, -s2])
-            assert allclose(sum(x.normal**2), 1)
+            assert numpy.allclose(x.normal, [-s2, -s2])
+            assert numpy.allclose(sum(x.normal**2), 1)
 
             msg = 'Triangle %d' %x.triangle_id + ' is not in %s' %(intersected_triangles)
@@ -1598 +1582 @@
         assert len(L) == 1
         assert L[0].triangle_id == 3
-        assert allclose(L[0].length, 0.5)        
-        assert allclose(L[0].normal, [-1,0])                
+        assert numpy.allclose(L[0].length, 0.5)        
+        assert numpy.allclose(L[0].normal, [-1,0])                
 
 
@@ -1610 +1594 @@
         assert len(L) == 1
         assert L[0].triangle_id == 3
-        assert allclose(L[0].length, 0.4)
-        assert allclose(L[0].normal, [-1,0])
+        assert numpy.allclose(L[0].length, 0.4)
+        assert numpy.allclose(L[0].normal, [-1,0])
 
         intersected_triangles = [3]
@@ -1623 +1607 @@
         assert len(L) == 1
         assert L[0].triangle_id == 3
-        assert allclose(L[0].length, 0.4)
-        assert allclose(L[0].normal, [1,0])                
+        assert numpy.allclose(L[0].length, 0.4)
+        assert numpy.allclose(L[0].normal, [1,0])                
         
 
@@ -1658 +1642 @@
         for x in L:
             if x.triangle_id == 3:
-                assert allclose(x.length, 0.5)        
-                assert allclose(x.normal, [-1,0])
+                assert numpy.allclose(x.length, 0.5)        
+                assert numpy.allclose(x.normal, [-1,0])
                 
             if x.triangle_id == 2:
-                assert allclose(x.length, s2)
-                assert allclose(x.normal, [-s2,-s2])
+                assert numpy.allclose(x.length, s2)
+                assert numpy.allclose(x.normal, [-s2,-s2])
 
 
@@ -1696 +1680 @@
         for x in L:
             if x.triangle_id == 3:
-                assert allclose(x.length, 0.5)        
-                assert allclose(x.normal, [-1,0])
+                assert numpy.allclose(x.length, 0.5)        
+                assert numpy.allclose(x.normal, [-1,0])
                 
             if x.triangle_id == 2:
                 msg = str(x.length)
-                assert allclose(x.length, s2), msg
-                assert allclose(x.normal, [-s2,-s2])
+                assert numpy.allclose(x.length, s2), msg
+                assert numpy.allclose(x.normal, [-s2,-s2])
 
             if x.triangle_id == 5:
-                segvec = array([line[2][0]-1,
-                                line[2][1]-1])
+                segvec = numpy.array([line[2][0]-1,
+                                      line[2][1]-1])
                 msg = str(x.length)
-                assert allclose(x.length, sqrt(sum(segvec**2))), msg
-                assert allclose(x.normal, [-s2,-s2])                                                
+                assert numpy.allclose(x.length, sqrt(sum(segvec**2))), msg
+                assert numpy.allclose(x.normal, [-s2,-s2])                                                
 
 
@@ -1747 +1731 @@
         for x in L:
             if x.triangle_id == 3:
-                assert allclose(x.length, 0.5)        
-                assert allclose(x.normal, [-1,0])
+                assert numpy.allclose(x.length, 0.5)        
+                assert numpy.allclose(x.normal, [-1,0])
                 
             if x.triangle_id == 2:
                 msg = str(x.length)
-                assert allclose(x.length, s2), msg
-                assert allclose(x.normal, [-s2,-s2])
+                assert numpy.allclose(x.length, s2), msg
+                assert numpy.allclose(x.normal, [-s2,-s2])
 
             if x.triangle_id == 5:
                 if x.segment == ((1.0, 1.0), (1.25, 0.75)):                    
-                    segvec = array([line[2][0]-1,
-                                    line[2][1]-1])
+                    segvec = numpy.array([line[2][0]-1,
+                                          line[2][1]-1])
                     msg = str(x.length)
-                    assert allclose(x.length, sqrt(sum(segvec**2))), msg
-                    assert allclose(x.normal, [-s2,-s2])
+                    assert numpy.allclose(x.length, sqrt(sum(segvec**2))), msg
+                    assert numpy.allclose(x.normal, [-s2,-s2])
                 elif x.segment == ((1.25, 0.75), (1.5, 1.0)):
-                    segvec = array([1.5-line[2][0],
-                                    1.0-line[2][1]])
+                    segvec = numpy.array([1.5-line[2][0],
+                                          1.0-line[2][1]])
                     
-                    assert allclose(x.length, sqrt(sum(segvec**2))), msg
-                    assert allclose(x.normal, [s2,-s2])
+                    assert numpy.allclose(x.length, sqrt(sum(segvec**2))), msg
+                    assert numpy.allclose(x.normal, [s2,-s2])
                 else:
                     msg = 'Unknow segment: %s' %x.segment
@@ -1775 +1759 @@
                     
             if x.triangle_id == 6:
-                assert allclose(x.normal, [s2,-s2])
-                assert allclose(x.segment, ((1.5, 1.0), (2, 1.5)))
+                assert numpy.allclose(x.normal, [s2,-s2])
+                assert numpy.allclose(x.segment, ((1.5, 1.0), (2, 1.5)))
 
 
@@ -1842 +1826 @@
             ref_length = line[1][1] - line[0][1]
             #print ref_length, total_length
-            assert allclose(total_length, ref_length)
+            assert numpy.allclose(total_length, ref_length)
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py

@@ -3 +3 @@
 
 import unittest
-
-from Numeric import allclose, array
-
+import numpy
 
 #from anuga.pyvolution.pmesh2domain import *
@@ -66 +64 @@
 
          tags = {}
-         b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
-         b2 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
-         b3 =  Dirichlet_boundary(conserved_quantities = array([2.0]))
+         b1 =  Dirichlet_boundary(conserved_quantities = numpy.array([0.0]))
+         b2 =  Dirichlet_boundary(conserved_quantities = numpy.array([1.0]))
+         b3 =  Dirichlet_boundary(conserved_quantities = numpy.array([2.0]))
          tags["1"] = b1
          tags["2"] = b2
@@ -80 +78 @@
          answer = [[0., 8., 0.],
                    [0., 10., 8.]]
-         assert allclose(domain.quantities['elevation'].vertex_values,
-                        answer)
+         assert numpy.allclose(domain.quantities['elevation'].vertex_values,
+                               answer)
 
          #print domain.quantities['stage'].vertex_values
          answer = [[0., 12., 10.],
                    [0., 10., 12.]]
-         assert allclose(domain.quantities['stage'].vertex_values,
-                        answer)
+         assert numpy.allclose(domain.quantities['stage'].vertex_values,
+                               answer)
 
          #print domain.quantities['friction'].vertex_values
          answer = [[0.01, 0.04, 0.03],
                    [0.01, 0.02, 0.04]]
-         assert allclose(domain.quantities['friction'].vertex_values,
-                        answer)
-
-         #print domain.quantities['friction'].vertex_values
-         assert allclose(domain.tagged_elements['dsg'][0],0)
-         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
+         assert numpy.allclose(domain.quantities['friction'].vertex_values,
+                               answer)
+
+         #print domain.quantities['friction'].vertex_values
+         assert numpy.allclose(domain.tagged_elements['dsg'][0],0)
+         assert numpy.allclose(domain.tagged_elements['ole nielsen'][0],1)
 
          self.failUnless( domain.boundary[(1, 0)]  == '1',
@@ -159 +157 @@
         
          tags = {}
-         b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
-         b2 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
-         b3 =  Dirichlet_boundary(conserved_quantities = array([2.0]))
+         b1 =  Dirichlet_boundary(conserved_quantities = numpy.array([0.0]))
+         b2 =  Dirichlet_boundary(conserved_quantities = numpy.array([1.0]))
+         b3 =  Dirichlet_boundary(conserved_quantities = numpy.array([2.0]))
          tags["1"] = b1
          tags["2"] = b2
@@ -174 +172 @@
          answer = [[0., 8., 0.],
                    [0., 10., 8.]]
-         assert allclose(domain.quantities['elevation'].vertex_values,
-                        answer)
+         assert numpy.allclose(domain.quantities['elevation'].vertex_values,
+                               answer)
 
          #print domain.quantities['stage'].vertex_values
          answer = [[0., 12., 10.],
                    [0., 10., 12.]]
-         assert allclose(domain.quantities['stage'].vertex_values,
-                        answer)
+         assert numpy.allclose(domain.quantities['stage'].vertex_values,
+                               answer)
 
          #print domain.quantities['friction'].vertex_values
          answer = [[0.01, 0.04, 0.03],
                    [0.01, 0.02, 0.04]]
-         assert allclose(domain.quantities['friction'].vertex_values,
-                        answer)
-
-         #print domain.quantities['friction'].vertex_values
-         assert allclose(domain.tagged_elements['dsg'][0],0)
-         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
+         assert numpy.allclose(domain.quantities['friction'].vertex_values,
+                               answer)
+
+         #print domain.quantities['friction'].vertex_values
+         assert numpy.allclose(domain.tagged_elements['dsg'][0],0)
+         assert numpy.allclose(domain.tagged_elements['ole nielsen'][0],1)
 
          self.failUnless( domain.boundary[(1, 0)]  == '1',
@@ -228 +226 @@
          #domain.set_tag_dict(tag_dict)
          #Boundary tests
-         b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
-         b2 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
-         b3 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
+         b1 =  Dirichlet_boundary(conserved_quantities = numpy.array([0.0]))
+         b2 =  Dirichlet_boundary(conserved_quantities = numpy.array([1.0]))
+         b3 =  Dirichlet_boundary(conserved_quantities = numpy.array([1.0]))
          #test adding a boundary
          tags = {}

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -7 +7 @@
 from quantity import *
 from anuga.config import epsilon
-from Numeric import allclose, array, ones, Float
+import numpy
 
 from anuga.fit_interpolate.fit import fit_to_mesh
@@ -18 +18 @@
 #Aux for fit_interpolate.fit example
 def linear_function(point):
-    point = array(point)
+    point = numpy.array(point)
     return point[:,0]+point[:,1]
 
@@ -63 +63 @@
 
         quantity = Quantity(self.mesh1, [[1,2,3]])
-        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
+        assert numpy.allclose(quantity.vertex_values, [[1.,2.,3.]])
 
         try:
@@ -84 +84 @@
 
         quantity = Quantity(self.mesh1)
-        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
-
-
-        quantity = Quantity(self.mesh4)
-        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
+        assert numpy.allclose(quantity.vertex_values, [[0.,0.,0.]])
+
+
+        quantity = Quantity(self.mesh4)
+        assert numpy.allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
                                                  [0.,0.,0.], [0.,0.,0.]])
 
@@ -94 +94 @@
     def test_interpolation(self):
         quantity = Quantity(self.mesh1, [[1,2,3]])
-        assert allclose(quantity.centroid_values, [2.0]) #Centroid
-
-        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
+        assert numpy.allclose(quantity.centroid_values, [2.0]) #Centroid
+
+        assert numpy.allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
 
 
@@ -102 +102 @@
         quantity = Quantity(self.mesh4,
                             [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
+        assert numpy.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
 
 
@@ -108 +108 @@
 
         #print quantity.vertex_values
-        assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],
-                                                 [3.+2./3, 6.+2./3, 4.+2./3],
-                                                 [4.6, 3.4, 1.],
-                                                 [-5.0, 1.0, 4.0]])
+        assert numpy.allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],
+                                                       [3.+2./3, 6.+2./3, 4.+2./3],
+                                                       [4.6, 3.4, 1.],
+                                                       [-5.0, 1.0, 4.0]])
 
         #print quantity.edge_values
-        assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25],
-                                               [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],
-                                               [2.2, 2.8, 4.0],
-                                               [2.5, -0.5, -2.0]])
-
-
-        #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
-        #                                       [5., 5., 5.],
-        #                                       [4.5, 4.5, 0.],
-        #                                       [3.0, -1.5, -1.5]])
+        assert numpy.allclose(quantity.edge_values, [[1.25, 3.5, 1.25],
+                                                     [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],
+                                                     [2.2, 2.8, 4.0],
+                                                     [2.5, -0.5, -2.0]])
+
+
+        #assert numpy.allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
+        #                                             [5., 5., 5.],
+        #                                             [4.5, 4.5, 0.],
+        #                                             [3.0, -1.5, -1.5]])
 
     def test_get_extrema_1(self):
         quantity = Quantity(self.mesh4,
                                       [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
+        assert numpy.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
 
         v = quantity.get_maximum_value()
@@ -148 +148 @@
 
         v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 5)
+        assert numpy.allclose(v, 5)
 
 
         x,y = quantity.get_minimum_location()
         v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 0)
+        assert numpy.allclose(v, 0)
 
 
@@ -192 +192 @@
 
         v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 6)
+        assert numpy.allclose(v, 6)
 
         x,y = quantity.get_minimum_location()        
         v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 2)
+        assert numpy.allclose(v, 2)
 
         #Multiple locations for maximum -
         #Test that the algorithm picks the first occurrence        
         v = quantity.get_maximum_value(indices=[0,1,2])
-        assert allclose(v, 4)
+        assert numpy.allclose(v, 4)
 
         i = quantity.get_maximum_index(indices=[0,1,2])
@@ -212 +212 @@
 
         v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 4)        
+        assert numpy.allclose(v, 4)        
 
         # More test of indices......
         v = quantity.get_maximum_value(indices=[2,3])
-        assert allclose(v, 6)
+        assert numpy.allclose(v, 6)
 
         i = quantity.get_maximum_index(indices=[2,3])
@@ -227 +227 @@
 
         v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 6)        
+        assert numpy.allclose(v, 6)        
 
         
@@ -244 +244 @@
         quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
                             location = 'vertices')
-        assert allclose(quantity.vertex_values,
-                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
-        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
-                                               [5., 5., 5.],
-                                               [4.5, 4.5, 0.],
-                                               [3.0, -1.5, -1.5]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+        assert numpy.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
+        assert numpy.allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
+                                                     [5., 5., 5.],
+                                                     [4.5, 4.5, 0.],
+                                                     [3.0, -1.5, -1.5]])
 
 
         # Test default
         quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-        assert allclose(quantity.vertex_values,
-                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
-        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
-                                               [5., 5., 5.],
-                                               [4.5, 4.5, 0.],
-                                               [3.0, -1.5, -1.5]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+        assert numpy.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
+        assert numpy.allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
+                                                     [5., 5., 5.],
+                                                     [4.5, 4.5, 0.],
+                                                     [3.0, -1.5, -1.5]])
 
         # Test centroids
         quantity.set_values([1,2,3,4], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
+        assert numpy.allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
 
         # Test exceptions
@@ -288 +288 @@
 
         quantity.set_values(1.0, location = 'vertices')
-        assert allclose(quantity.vertex_values,
-                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
-
-        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
-        assert allclose(quantity.edge_values, [[1, 1, 1],
-                                               [1, 1, 1],
-                                               [1, 1, 1],
-                                               [1, 1, 1]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
+
+        assert numpy.allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
+        assert numpy.allclose(quantity.edge_values, [[1, 1, 1],
+                                                     [1, 1, 1],
+                                                     [1, 1, 1],
+                                                     [1, 1, 1]])
 
 
         quantity.set_values(2.0, location = 'centroids')
-        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
+        assert numpy.allclose(quantity.centroid_values, [2, 2, 2, 2])
 
 
@@ -310 +310 @@
         quantity.set_values(f, location = 'vertices')
         #print "quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values,
-                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
-        assert allclose(quantity.centroid_values,
-                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
-        assert allclose(quantity.edge_values,
-                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
+        assert numpy.allclose(quantity.centroid_values,
+                              [4.0/3, 8.0/3, 10.0/3, 10.0/3])
+        assert numpy.allclose(quantity.edge_values,
+                              [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
 
 
         quantity.set_values(f, location = 'centroids')
-        assert allclose(quantity.centroid_values,
-                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
+        assert numpy.allclose(quantity.centroid_values,
+                              [4.0/3, 8.0/3, 10.0/3, 10.0/3])
 
 
@@ -331 +331 @@
         #print 'Q', quantity.get_integral()
 
-        assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
+        assert numpy.allclose(quantity.get_integral(), self.mesh4.get_area() * const)
 
         #Try with a linear function
@@ -342 +342 @@
         ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
 
-        assert allclose (quantity.get_integral(), ref_integral)
+        assert numpy.allclose (quantity.get_integral(), ref_integral)
 
 
@@ -350 +350 @@
         quantity.set_vertex_values([0,1,2,3,4,5])
 
-        assert allclose(quantity.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
-        assert allclose(quantity.centroid_values,
-                        [1., 7./3, 11./3, 8./3]) #Centroid
-        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
-                                               [3., 2.5, 1.5],
-                                               [3.5, 4.5, 3.],
-                                               [2.5, 3.5, 2]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        assert numpy.allclose(quantity.centroid_values,
+                              [1., 7./3, 11./3, 8./3]) #Centroid
+        assert numpy.allclose(quantity.edge_values, [[1., 1.5, 0.5],
+                                                     [3., 2.5, 1.5],
+                                                     [3.5, 4.5, 3.],
+                                                     [2.5, 3.5, 2]])
 
 
@@ -365 +365 @@
         quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
 
-        assert allclose(quantity.vertex_values,
-                        [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
 
 
@@ -376 +376 @@
 
 
-        assert allclose(quantity.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
         #Centroid
-        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
-
-        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
-                                               [3., 2.5, 1.5],
-                                               [3.5, 4.5, 3.],
-                                               [2.5, 3.5, 2]])
+        assert numpy.allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
+
+        assert numpy.allclose(quantity.edge_values, [[1., 1.5, 0.5],
+                                                     [3., 2.5, 1.5],
+                                                     [3.5, 4.5, 3.],
+                                                     [2.5, 3.5, 2]])
 
 
@@ -399 +399 @@
 
         quantity.set_values([0,2,3,5], indices=[0,2,3,5])
-        assert allclose(quantity.vertex_values,
-                        [[0,0,2], [0,2,0], [0,2,5], [3,0,0]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[0,0,2], [0,2,0], [0,2,5], [3,0,0]])
 
 
@@ -408 +408 @@
 
         # Indices refer to triangle numbers
-        assert allclose(quantity.vertex_values,
-                        [[3.14,3.14,3.14], [0,0,0],
-                         [3.14,3.14,3.14], [0,0,0]])        
+        assert numpy.allclose(quantity.vertex_values,
+                              [[3.14,3.14,3.14], [0,0,0],
+                               [3.14,3.14,3.14], [0,0,0]])        
         
 
@@ -419 +419 @@
         quantity.set_values(3.14, polygon=polygon)
         
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0], [0,0,0], [0,0,0],
-                         [3.14,3.14,3.14]])                
+        assert numpy.allclose(quantity.vertex_values,
+                              [[0,0,0], [0,0,0], [0,0,0],
+                               [3.14,3.14,3.14]])                
 
 
@@ -429 +429 @@
         quantity.set_values(0.0)
         quantity.set_values(3.14, location='centroids', polygon=polygon)
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0],
-                         [3.14,3.14,3.14],
-                         [3.14,3.14,3.14],                         
-                         [0,0,0]])                
+        assert numpy.allclose(quantity.vertex_values,
+                              [[0,0,0],
+                               [3.14,3.14,3.14],
+                               [3.14,3.14,3.14],                         
+                               [0,0,0]])                
 
 
@@ -441 +441 @@
         #print 'Here 2'
         quantity.set_values(3.14, polygon=polygon)
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0],
-                         [3.14,3.14,3.14],
-                         [3.14,3.14,3.14],                         
-                         [0,0,0]])                
+        assert numpy.allclose(quantity.vertex_values,
+                              [[0,0,0],
+                               [3.14,3.14,3.14],
+                               [3.14,3.14,3.14],                         
+                               [0,0,0]])                
         
 
@@ -478 +478 @@
 
         # Indices refer to triangle numbers here - not vertices (why?)
-        assert allclose(quantity.vertex_values,
-                        [[3.14,3.14,3.14], [0,0,0],
-                         [3.14,3.14,3.14], [0,0,0]])        
+        assert numpy.allclose(quantity.vertex_values,
+                              [[3.14,3.14,3.14], [0,0,0],
+                               [3.14,3.14,3.14], [0,0,0]])        
         
 
 
         # Now try with polygon (pick points where y>2)
-        polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]])
+        polygon = numpy.array([[0,2.1], [4,2.1], [4,7], [0,7]])
         polygon += [G.xllcorner, G.yllcorner]
         
@@ -491 +491 @@
         quantity.set_values(3.14, polygon=polygon, location='centroids')
         
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0], [0,0,0], [0,0,0],
-                         [3.14,3.14,3.14]])                
+        assert numpy.allclose(quantity.vertex_values,
+                              [[0,0,0], [0,0,0], [0,0,0],
+                               [3.14,3.14,3.14]])                
 
 
         # Another polygon (pick triangle 1 and 2 (rightmost triangles)
-        polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])
+        polygon = numpy.array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])
         polygon += [G.xllcorner, G.yllcorner]
         
@@ -503 +503 @@
         quantity.set_values(3.14, polygon=polygon)
 
-        assert allclose(quantity.vertex_values,
-                        [[0,0,0],
-                         [3.14,3.14,3.14],
-                         [3.14,3.14,3.14],                         
-                         [0,0,0]])                
+        assert numpy.allclose(quantity.vertex_values,
+                              [[0,0,0],
+                               [3.14,3.14,3.14],
+                               [3.14,3.14,3.14],                         
+                               [0,0,0]])                
 
 
@@ -540 +540 @@
         answer = linear_function(quantity.domain.get_vertex_coordinates())
         #print quantity.vertex_values, answer
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
@@ -553 +553 @@
         #print vertex_attributes
         quantity.set_values(vertex_attributes)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
@@ -593 +593 @@
                           alpha = 0)
 
-        assert allclose( ref, [0,5,5] )
+        assert numpy.allclose( ref, [0,5,5] )
 
 
@@ -610 +610 @@
         #                    data_georef = data_georef,
         #                    alpha = 0)
-        assert allclose(quantity.vertex_values.flat, ref)
+        assert numpy.allclose(quantity.vertex_values.ravel(), ref)
 
 
@@ -621 +621 @@
 
         quantity.set_values(geospatial_data = geo, alpha = 0)
-        assert allclose(quantity.vertex_values.flat, ref)
+        assert numpy.allclose(quantity.vertex_values.ravel(), ref)
 
 
@@ -666 +666 @@
         answer = linear_function(quantity.domain.get_vertex_coordinates())
 
-        #print quantity.vertex_values.flat
+        #print quantity.vertex_values.ravel()
         #print answer
 
 
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
         #Check that values can be set from file using default attribute
         quantity.set_values(filename = ptsfile, alpha = 0)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         #Cleanup
@@ -730 +730 @@
         #print self.mesh4.nodes
         #print inside_polygon(self.mesh4.nodes, polygon)
-        assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)
+        assert numpy.allclose(inside_polygon(self.mesh4.nodes, polygon), 4)
 
         #print quantity.domain.get_vertex_coordinates()
@@ -748 +748 @@
         answer = linear_function(points)
 
-        #print quantity.vertex_values.flat
+        #print quantity.vertex_values.ravel()
         #print answer
 
         # Check vertices in polygon have been set
-        assert allclose(take(quantity.vertex_values.flat, indices),
-                             answer)
+        assert numpy.allclose(take(quantity.vertex_values.ravel(), indices),
+                                   answer)
 
         # Check vertices outside polygon are zero
         indices = outside_polygon(quantity.domain.get_vertex_coordinates(), 
                                   polygon)        
-        assert allclose(take(quantity.vertex_values.flat, indices),
-                             0.0)        
+        assert numpy.allclose(take(quantity.vertex_values.ravel(), indices),
+                                   0.0)        
 
         #Cleanup
@@ -815 +815 @@
                             verbose=False)
         answer = linear_function(quantity.domain.get_vertex_coordinates())
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
@@ -821 +821 @@
         quantity.set_values(filename=ptsfile, 
                             alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         # Check cache
@@ -866 +866 @@
         answer = linear_function(quantity.domain.get_vertex_coordinates())
 
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
+        #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel()
         #print "answer",answer
 
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
         #Check that values can be set from file using default attribute
         quantity.set_values(filename=txt_file, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         #Cleanup
@@ -910 +910 @@
         answer = linear_function(quantity.domain.get_vertex_coordinates())
 
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
+        #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel()
         #print "answer",answer
 
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
         #Check that values can be set from file using default attribute
         quantity.set_values(filename=txt_file, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         #Cleanup
@@ -957 +957 @@
                                       'vertices', None)
         answer = linear_function(quantity.domain.get_vertex_coordinates())
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
+        #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel()
         #print "answer",answer
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         #Check that values can be set from file
@@ -965 +965 @@
                             attribute_name=att, alpha=0)
         answer = linear_function(quantity.domain.get_vertex_coordinates())
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
+        #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel()
         #print "answer",answer
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
         #Check that values can be set from file using default attribute
         quantity.set_values(filename=txt_file, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         #Cleanup
@@ -1039 +1039 @@
                                       max_read_lines=2)
         answer = linear_function(quantity.domain.get_vertex_coordinates())
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
+        #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel()
         #print "answer",answer
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         #Check that values can be set from file
@@ -1047 +1047 @@
                             attribute_name=att, alpha=0)
         answer = linear_function(quantity.domain.get_vertex_coordinates())
-        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
+        #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel()
         #print "answer",answer
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
         #Check that values can be set from file using default attribute
         quantity.set_values(filename=txt_file, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         #Cleanup
@@ -1128 +1128 @@
         answer = linear_function(quantity.domain.get_vertex_coordinates())
 
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
         #Check that values can be set from file using default attribute
         quantity.set_values(filename=ptsfile, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         #Cleanup
@@ -1207 +1207 @@
 
 
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
         #Check that values can be set from file using default attribute
         quantity.set_values(filename=ptsfile, alpha=0)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
         #Cleanup
@@ -1226 +1226 @@
         quantity1.set_vertex_values([0,1,2,3,4,5])
 
-        assert allclose(quantity1.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        assert numpy.allclose(quantity1.vertex_values,
+                              [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
 
         quantity2 = Quantity(self.mesh4)
         quantity2.set_values(quantity=quantity1)
-        assert allclose(quantity2.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        assert numpy.allclose(quantity2.vertex_values,
+                              [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
         quantity2.set_values(quantity = 2*quantity1)
-        assert allclose(quantity2.vertex_values,
-                        [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
+        assert numpy.allclose(quantity2.vertex_values,
+                              [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
 
         quantity2.set_values(quantity = 2*quantity1 + 3)
-        assert allclose(quantity2.vertex_values,
-                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
+        assert numpy.allclose(quantity2.vertex_values,
+                              [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
 
 
         #Check detection of quantity as first orgument
         quantity2.set_values(2*quantity1 + 3)
-        assert allclose(quantity2.vertex_values,
-                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
+        assert numpy.allclose(quantity2.vertex_values,
+                              [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
 
 
@@ -1262 +1262 @@
         polygon = [[1.0, 1.0], [4.0, 1.0],
                    [4.0, 4.0], [1.0, 4.0]]
-        assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)                   
+        assert numpy.allclose(inside_polygon(self.mesh4.nodes, polygon), 4)                   
         
         quantity1 = Quantity(self.mesh4)
         quantity1.set_vertex_values([0,1,2,3,4,5])
 
-        assert allclose(quantity1.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        assert numpy.allclose(quantity1.vertex_values,
+                              [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
 
@@ -1276 +1276 @@
                              
         msg = 'Only node #4(e) at (2,2) should have values applied '
-        assert allclose(quantity2.vertex_values,
-                        [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg        
-                        #bac,     bce,     ecf,     dbe
+        assert numpy.allclose(quantity2.vertex_values,
+                              [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg        
+                              #bac,     bce,     ecf,     dbe
                         
 
@@ -1287 +1287 @@
         quantity1.set_vertex_values([0,1,2,3,4,5])
 
-        assert allclose(quantity1.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        assert numpy.allclose(quantity1.vertex_values,
+                              [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
 
@@ -1304 +1304 @@
         # Negation
         Q = -quantity1
-        assert allclose(Q.vertex_values, -quantity1.vertex_values)
-        assert allclose(Q.centroid_values, -quantity1.centroid_values)
-        assert allclose(Q.edge_values, -quantity1.edge_values)
+        assert numpy.allclose(Q.vertex_values, -quantity1.vertex_values)
+        assert numpy.allclose(Q.centroid_values, -quantity1.centroid_values)
+        assert numpy.allclose(Q.edge_values, -quantity1.edge_values)
 
         # Addition
         Q = quantity1 + 7
-        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
-        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
-        assert allclose(Q.edge_values, quantity1.edge_values + 7)
+        assert numpy.allclose(Q.vertex_values, quantity1.vertex_values + 7)
+        assert numpy.allclose(Q.centroid_values, quantity1.centroid_values + 7)
+        assert numpy.allclose(Q.edge_values, quantity1.edge_values + 7)
 
         Q = 7 + quantity1
-        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
-        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
-        assert allclose(Q.edge_values, quantity1.edge_values + 7)
+        assert numpy.allclose(Q.vertex_values, quantity1.vertex_values + 7)
+        assert numpy.allclose(Q.centroid_values, quantity1.centroid_values + 7)
+        assert numpy.allclose(Q.edge_values, quantity1.edge_values + 7)
 
         Q = quantity1 + quantity2
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values + quantity2.vertex_values)
-        assert allclose(Q.centroid_values,
-                        quantity1.centroid_values + quantity2.centroid_values)
-        assert allclose(Q.edge_values,
-                        quantity1.edge_values + quantity2.edge_values)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values + quantity2.vertex_values)
+        assert numpy.allclose(Q.centroid_values,
+                              quantity1.centroid_values + quantity2.centroid_values)
+        assert numpy.allclose(Q.edge_values,
+                              quantity1.edge_values + quantity2.edge_values)
 
 
         Q = quantity1 + quantity2 - 3
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values + quantity2.vertex_values - 3)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values + quantity2.vertex_values - 3)
 
         Q = quantity1 - quantity2
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values - quantity2.vertex_values)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values - quantity2.vertex_values)
 
         #Scaling
         Q = quantity1*3
-        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
-        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
-        assert allclose(Q.edge_values, quantity1.edge_values*3)
+        assert numpy.allclose(Q.vertex_values, quantity1.vertex_values*3)
+        assert numpy.allclose(Q.centroid_values, quantity1.centroid_values*3)
+        assert numpy.allclose(Q.edge_values, quantity1.edge_values*3)
         Q = 3*quantity1
-        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
+        assert numpy.allclose(Q.vertex_values, quantity1.vertex_values*3)
 
         #Multiplication
@@ -1351 +1351 @@
         #print quantity2.centroid_values
 
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values * quantity2.vertex_values)
 
         #Linear combinations
         Q = 4*quantity1 + 2
-        assert allclose(Q.vertex_values,
-                        4*quantity1.vertex_values + 2)
+        assert numpy.allclose(Q.vertex_values,
+                              4*quantity1.vertex_values + 2)
 
         Q = quantity1*quantity2 + 2
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values + 2)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values * quantity2.vertex_values + 2)
 
         Q = quantity1*quantity2 + quantity3
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values +
-                        quantity3.vertex_values)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values * quantity2.vertex_values +
+                              quantity3.vertex_values)
         Q = quantity1*quantity2 + 3*quantity3
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values +
-                        3*quantity3.vertex_values)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values * quantity2.vertex_values +
+                             3*quantity3.vertex_values)
         Q = quantity1*quantity2 + 3*quantity3 + 5.0
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values +
-                        3*quantity3.vertex_values + 5)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values * quantity2.vertex_values +
+                              3*quantity3.vertex_values + 5)
 
         Q = quantity1*quantity2 - quantity3
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values * quantity2.vertex_values -
-                        quantity3.vertex_values)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values * quantity2.vertex_values -
+                              quantity3.vertex_values)
         Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
-        assert allclose(Q.vertex_values,
-                        1.5*quantity1.vertex_values * quantity2.vertex_values -
-                        3*quantity3.vertex_values + 5)
+        assert numpy.allclose(Q.vertex_values,
+                              1.5*quantity1.vertex_values * quantity2.vertex_values -
+                              3*quantity3.vertex_values + 5)
 
         #Try combining quantities and arrays and scalars
-        Q = 1.5*quantity1*quantity2.vertex_values -\
+        Q = 1.5*quantity1*quantity2.vertex_values - \
             3*quantity3.vertex_values + 5.0
-        assert allclose(Q.vertex_values,
-                        1.5*quantity1.vertex_values * quantity2.vertex_values -
-                        3*quantity3.vertex_values + 5)
+        assert numpy.allclose(Q.vertex_values,
+                              1.5*quantity1.vertex_values * quantity2.vertex_values -
+                              3*quantity3.vertex_values + 5)
 
 
         #Powers
         Q = quantity1**2
-        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
+        assert numpy.allclose(Q.vertex_values, quantity1.vertex_values**2)
 
         Q = quantity1**2 +quantity2**2
-        assert allclose(Q.vertex_values,
-                        quantity1.vertex_values**2 + \
-                        quantity2.vertex_values**2)
+        assert numpy.allclose(Q.vertex_values,
+                              quantity1.vertex_values**2 +
+                              quantity2.vertex_values**2)
 
         Q = (quantity1**2 +quantity2**2)**0.5
-        assert allclose(Q.vertex_values,
-                        (quantity1.vertex_values**2 + \
-                        quantity2.vertex_values**2)**0.5)
+        assert numpy.allclose(Q.vertex_values,
+                              (quantity1.vertex_values**2 + 
+                               quantity2.vertex_values**2)**0.5)
 
 
@@ -1431 +1431 @@
         #The central triangle (1)
         #(using standard gradient based on neigbours controid values)
-        assert allclose(a[1], 2.0)
-        assert allclose(b[1], 0.0)
+        assert numpy.allclose(a[1], 2.0)
+        assert numpy.allclose(b[1], 0.0)
 
 
@@ -1438 +1438 @@
         #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
         #2  = 4  + a*(-2/3)  + b*(-2/3)
-        assert allclose(a[0] + b[0], 3)
+        assert numpy.allclose(a[0] + b[0], 3)
         #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
-        assert allclose(a[0] - b[0], 0)
+        assert numpy.allclose(a[0] - b[0], 0)
 
 
@@ -1446 +1446 @@
         #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
         #6  = 4  + a*(4/3)  + b*(-2/3)
-        assert allclose(2*a[2] - b[2], 3)
+        assert numpy.allclose(2*a[2] - b[2], 3)
         #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
-        assert allclose(a[2] + 2*b[2], 0)
+        assert numpy.allclose(a[2] + 2*b[2], 0)
 
 
@@ -1454 +1454 @@
         #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
         #2  = 4  + a*(-2/3)  + b*(4/3)
-        assert allclose(a[3] - 2*b[3], 3)
+        assert numpy.allclose(a[3] - 2*b[3], 3)
         #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
-        assert allclose(2*a[3] + b[3], 0)
+        assert numpy.allclose(2*a[3] + b[3], 0)
 
 
@@ -1464 +1464 @@
 
         #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
-        assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
-        assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
+        assert numpy.allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
+        assert numpy.allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
 
 
         #a = 1.2, b=-0.6
         #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
-        assert allclose(quantity.vertex_values[2,2], 8)
+        assert numpy.allclose(quantity.vertex_values[2,2], 8)
 
     def test_get_gradients(self):
@@ -1489 +1489 @@
         #The central triangle (1)
         #(using standard gradient based on neigbours controid values)
-        assert allclose(a[1], 2.0)
-        assert allclose(b[1], 0.0)
+        assert numpy.allclose(a[1], 2.0)
+        assert numpy.allclose(b[1], 0.0)
 
 
@@ -1496 +1496 @@
         #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
         #2  = 4  + a*(-2/3)  + b*(-2/3)
-        assert allclose(a[0] + b[0], 3)
+        assert numpy.allclose(a[0] + b[0], 3)
         #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
-        assert allclose(a[0] - b[0], 0)
+        assert numpy.allclose(a[0] - b[0], 0)
 
 
@@ -1504 +1504 @@
         #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
         #6  = 4  + a*(4/3)  + b*(-2/3)
-        assert allclose(2*a[2] - b[2], 3)
+        assert numpy.allclose(2*a[2] - b[2], 3)
         #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
-        assert allclose(a[2] + 2*b[2], 0)
+        assert numpy.allclose(a[2] + 2*b[2], 0)
 
 
@@ -1512 +1512 @@
         #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
         #2  = 4  + a*(-2/3)  + b*(4/3)
-        assert allclose(a[3] - 2*b[3], 3)
+        assert numpy.allclose(a[3] - 2*b[3], 3)
         #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
-        assert allclose(2*a[3] + b[3], 0)
+        assert numpy.allclose(2*a[3] + b[3], 0)
 
 
@@ -1532 +1532 @@
         #print a, b
 
-        assert allclose(a[1], 3.0)
-        assert allclose(b[1], 1.0)
+        assert numpy.allclose(a[1], 3.0)
+        assert numpy.allclose(b[1], 1.0)
 
         #Work out the others
@@ -1540 +1540 @@
 
         #print quantity.vertex_values
-        assert allclose(quantity.vertex_values[1,0], 2.0)
-        assert allclose(quantity.vertex_values[1,1], 6.0)
-        assert allclose(quantity.vertex_values[1,2], 8.0)
+        assert numpy.allclose(quantity.vertex_values[1,0], 2.0)
+        assert numpy.allclose(quantity.vertex_values[1,1], 6.0)
+        assert numpy.allclose(quantity.vertex_values[1,2], 8.0)
 
 
@@ -1550 +1550 @@
 
         #Set up for a gradient of (3,1), f(x) = 3x+y
-        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
-        d_values = array([1.0, 2.0, 3.0, 4.0])
+        c_values = numpy.array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
+        d_values = numpy.array([1.0, 2.0, 3.0, 4.0])
         quantity.set_values(c_values, location = 'centroids')
 
@@ -1558 +1558 @@
 
         #print quantity.vertex_values
-        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
+        assert numpy.allclose(quantity.centroid_values, quantity.centroid_backup_values)
 
 
@@ -1574 +1574 @@
         #Test centroids
         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+        assert numpy.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
 
         #Extrapolate
@@ -1581 +1581 @@
         #Check that gradient is zero
         a,b = quantity.get_gradients()
-        assert allclose(a, [0,0,0,0])
-        assert allclose(b, [0,0,0,0])
+        assert numpy.allclose(a, [0,0,0,0])
+        assert numpy.allclose(b, [0,0,0,0])
 
         #Check vertices but not edge values
-        assert allclose(quantity.vertex_values,
-                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
 
 
@@ -1614 +1614 @@
 
         #Assert that quantities are conserved
-        from Numeric import sum
         for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
+            assert numpy.allclose (quantity.centroid_values[k],
+                                   numpy.sum(quantity.vertex_values[k,:])/3)
 
 
@@ -1646 +1645 @@
 
         #Assert that quantities are conserved
-        from Numeric import sum
         for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
+            assert numpy.allclose (quantity.centroid_values[k],
+                                   numpy.sum(quantity.vertex_values[k,:])/3)
 
 
@@ -1676 +1674 @@
 
         #Assert that quantities are conserved
-        from Numeric import sum
         for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
+            assert numpy.allclose (quantity.centroid_values[k],
+                                   numpy.sum(quantity.vertex_values[k,:])/3)
 
 
@@ -1705 +1702 @@
 
         #Assert that quantities are conserved
-        from Numeric import sum
         for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
+            assert numpy.allclose (quantity.centroid_values[k],
+                                   numpy.sum(quantity.vertex_values[k,:])/3)
 
     def test_limiter2(self):
@@ -1718 +1714 @@
         #Test centroids
         quantity.set_values([2.,4.,8.,2.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
+        assert numpy.allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
 
 
@@ -1724 +1720 @@
         quantity.extrapolate_second_order()
 
-        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
+        assert numpy.allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
 
         #Limit
@@ -1731 +1727 @@
         # limited value for beta_w = 0.9
         
-        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
+        assert numpy.allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
         # limited values for beta_w = 0.5
-        #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
+        #assert numpy.allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
 
 
         #Assert that quantities are conserved
-        from Numeric import sum
         for k in range(quantity.centroid_values.shape[0]):
-            assert allclose (quantity.centroid_values[k],
-                             sum(quantity.vertex_values[k,:])/3)
+            assert numpy.allclose (quantity.centroid_values[k],
+                                   numpy.sum(quantity.vertex_values[k,:])/3)
 
 
@@ -1751 +1746 @@
         #Test centroids
         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+        assert numpy.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
 
 
@@ -1760 +1755 @@
         #quantity.interpolate_from_vertices_to_edges()
 
-        assert allclose(quantity.vertex_values,
-                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
-        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
-                                               [3,3,3], [4, 4, 4]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
+        assert numpy.allclose(quantity.edge_values, [[1,1,1], [2,2,2],
+                                                     [3,3,3], [4, 4, 4]])
 
 
@@ -1769 +1764 @@
         quantity = Quantity(self.mesh4)
 
-        quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)
+        quantity.vertex_values = numpy.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], numpy.float)
 
         quantity.interpolate_from_vertices_to_edges()
 
-        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
-                                               [3., 2.5, 1.5],
-                                               [3.5, 4.5, 3.],
-                                               [2.5, 3.5, 2]])
+        assert numpy.allclose(quantity.edge_values, [[1., 1.5, 0.5],
+                                                     [3., 2.5, 1.5],
+                                                     [3.5, 4.5, 3.],
+                                                     [2.5, 3.5, 2]])
 
 
@@ -1782 +1777 @@
         quantity = Quantity(self.mesh4)
 
-        quantity.edge_values = array([[1., 1.5, 0.5],
-                                [3., 2.5, 1.5],
-                                [3.5, 4.5, 3.],
-                                [2.5, 3.5, 2]],Float)
+        quantity.edge_values = numpy.array([[1., 1.5, 0.5],
+                                            [3., 2.5, 1.5],
+                                            [3.5, 4.5, 3.],
+                                            [2.5, 3.5, 2]], numpy.float)
 
         quantity.interpolate_from_edges_to_vertices()
 
-        assert allclose(quantity.vertex_values,
-                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        assert numpy.allclose(quantity.vertex_values,
+                              [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
 
 
@@ -1799 +1794 @@
         #Test centroids
         quantity.set_values([2.,4.,8.,2.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
+        assert numpy.allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
 
 
@@ -1805 +1800 @@
         quantity.extrapolate_second_order()
 
-        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
+        assert numpy.allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
 
 
@@ -1813 +1808 @@
         #Test centroids
         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+        assert numpy.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
 
         #Set explicit_update
-        quantity.explicit_update = array( [1.,1.,1.,1.] )
+        quantity.explicit_update = numpy.array( [1.,1.,1.,1.] )
 
         #Update with given timestep
         quantity.update(0.1)
 
-        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
-        assert allclose( quantity.centroid_values, x)
+        x = numpy.array([1, 2, 3, 4]) + numpy.array( [.1,.1,.1,.1] )
+        assert numpy.allclose( quantity.centroid_values, x)
 
     def test_update_semi_implicit(self):
@@ -1829 +1824 @@
         #Test centroids
         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+        assert numpy.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
 
         #Set semi implicit update
-        quantity.semi_implicit_update = array([1.,1.,1.,1.])
+        quantity.semi_implicit_update = numpy.array([1.,1.,1.,1.])
 
         #Update with given timestep
@@ -1838 +1833 @@
         quantity.update(timestep)
 
-        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
-        denom = ones(4, Float)-timestep*sem
-
-        x = array([1, 2, 3, 4])/denom
-        assert allclose( quantity.centroid_values, x)
+        sem = numpy.array([1.,1.,1.,1.])/numpy.array([1, 2, 3, 4])
+        denom = numpy.ones(4, numpy.float)-timestep*sem
+
+        x = numpy.array([1, 2, 3, 4])/denom
+        assert numpy.allclose( quantity.centroid_values, x)
 
 
@@ -1850 +1845 @@
         #Test centroids
         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+        assert numpy.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
 
         #Set explicit_update
-        quantity.explicit_update = array( [4.,3.,2.,1.] )
+        quantity.explicit_update = numpy.array( [4.,3.,2.,1.] )
 
         #Set semi implicit update
-        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
+        quantity.semi_implicit_update = numpy.array( [1.,1.,1.,1.] )
 
         #Update with given timestep
@@ -1862 +1857 @@
         quantity.update(0.1)
 
-        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
-        denom = ones(4, Float)-timestep*sem
-
-        x = array([1., 2., 3., 4.])
+        sem = numpy.array([1.,1.,1.,1.])/numpy.array([1, 2, 3, 4])
+        denom = numpy.ones(4, numpy.float)-timestep*sem
+
+        x = numpy.array([1., 2., 3., 4.])
         x /= denom
-        x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
-
-        assert allclose( quantity.centroid_values, x)
+        x += timestep*numpy.array( [4.0, 3.0, 2.0, 1.0] )
+
+        assert numpy.allclose( quantity.centroid_values, x)
 
 
@@ -1879 +1874 @@
         from mesh_factory import rectangular
         from shallow_water import Domain, Transmissive_boundary
-        from Numeric import zeros, Float
         from anuga.utilities.numerical_tools import mean
 
@@ -1906 +1900 @@
 
         bed = domain.quantities['elevation'].vertex_values
-        stage = zeros(bed.shape, Float)
+        stage = numpy.zeros(bed.shape, numpy.float)
 
         h = 0.03
@@ -1929 +1923 @@
 
         #First four points
-        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
-        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
-        assert allclose(A[2], Q[3,0])
-        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
+        assert numpy.allclose(A[0], (Q[0,2] + Q[1,1])/2)
+        assert numpy.allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
+        assert numpy.allclose(A[2], Q[3,0])
+        assert numpy.allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
 
         #Center point
-        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
-                               Q[5,0] + Q[6,2] + Q[7,1])/6)
+        assert numpy.allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
+                                     Q[5,0] + Q[6,2] + Q[7,1])/6)
 
 
         #Check V
-        assert allclose(V[0,:], [3,4,0])
-        assert allclose(V[1,:], [1,0,4])
-        assert allclose(V[2,:], [4,5,1])
-        assert allclose(V[3,:], [2,1,5])
-        assert allclose(V[4,:], [6,7,3])
-        assert allclose(V[5,:], [4,3,7])
-        assert allclose(V[6,:], [7,8,4])
-        assert allclose(V[7,:], [5,4,8])
+        assert numpy.allclose(V[0,:], [3,4,0])
+        assert numpy.allclose(V[1,:], [1,0,4])
+        assert numpy.allclose(V[2,:], [4,5,1])
+        assert numpy.allclose(V[3,:], [2,1,5])
+        assert numpy.allclose(V[4,:], [6,7,3])
+        assert numpy.allclose(V[5,:], [4,3,7])
+        assert numpy.allclose(V[6,:], [7,8,4])
+        assert numpy.allclose(V[7,:], [5,4,8])
 
         #Get smoothed stage with XY
         X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
 
-        assert allclose(A, A1)
-        assert allclose(V, V1)
+        assert numpy.allclose(A, A1)
+        assert numpy.allclose(V, V1)
 
         #Check XY
-        assert allclose(X[4], 0.5)
-        assert allclose(Y[4], 0.5)
-
-        assert allclose(X[7], 1.0)
-        assert allclose(Y[7], 0.5)
+        assert numpy.allclose(X[4], 0.5)
+        assert numpy.allclose(Y[4], 0.5)
+
+        assert numpy.allclose(X[7], 1.0)
+        assert numpy.allclose(Y[7], 0.5)
 
 
@@ -1969 +1963 @@
         from mesh_factory import rectangular
         from shallow_water import Domain, Transmissive_boundary
-        from Numeric import zeros, Float
         from anuga.utilities.numerical_tools import mean
 
@@ -1991 +1984 @@
 
         bed = domain.quantities['elevation'].vertex_values
-        stage = zeros(bed.shape, Float)
+        stage = numpy.zeros(bed.shape, numpy.float)
 
         h = 0.03
@@ -2005 +1998 @@
         stage = domain.quantities['stage']
         A, V = stage.get_vertex_values(xy=False, smooth=False)
-        Q = stage.vertex_values.flat
+        Q = stage.vertex_values.ravel()
 
         for k in range(8):
-            assert allclose(A[k], Q[k])
+            assert numpy.allclose(A[k], Q[k])
 
 
@@ -2021 +2014 @@
 
 
-        assert allclose(A, A1)
-        assert allclose(V, V1)
+        assert numpy.allclose(A, A1)
+        assert numpy.allclose(V, V1)
 
         #Check XY
-        assert allclose(X[1], 0.5)
-        assert allclose(Y[1], 0.5)
-        assert allclose(X[4], 0.0)
-        assert allclose(Y[4], 0.0)
-        assert allclose(X[12], 1.0)
-        assert allclose(Y[12], 0.0)
+        assert numpy.allclose(X[1], 0.5)
+        assert numpy.allclose(Y[1], 0.5)
+        assert numpy.allclose(X[4], 0.0)
+        assert numpy.allclose(Y[4], 0.0)
+        assert numpy.allclose(X[12], 1.0)
+        assert numpy.allclose(Y[12], 0.0)
 
 
@@ -2038 +2031 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2054 +2046 @@
         #print "quantity.centroid_values",quantity.centroid_values
 
-        assert allclose(quantity.centroid_values, [1,7])
+        assert numpy.allclose(quantity.centroid_values, [1,7])
 
         quantity.set_array_values([15,20,25], indices = indices)
-        assert allclose(quantity.centroid_values, [1,20])
+        assert numpy.allclose(quantity.centroid_values, [1,20])
 
         quantity.set_array_values([15,20,25], indices = indices)
-        assert allclose(quantity.centroid_values, [1,20])
+        assert numpy.allclose(quantity.centroid_values, [1,20])
 
     def test_setting_some_vertex_values(self):
@@ -2068 +2060 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2087 +2078 @@
                             indices = indices)
         #print "quantity.centroid_values",quantity.centroid_values
-        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
+        assert numpy.allclose(quantity.centroid_values, [1,7,3,4,5,6])
         
         value = [7]
@@ -2095 +2086 @@
                             indices = indices)
         #print "quantity.centroid_values",quantity.centroid_values
-        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
+        assert numpy.allclose(quantity.centroid_values, [1,7,3,4,5,6])
 
         value = [[15,20,25]]
         quantity.set_values(value, indices = indices)
         #print "1 quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values[1], value[0])
+        assert numpy.allclose(quantity.vertex_values[1], value[0])
 
 
@@ -2107 +2098 @@
         quantity.set_values(values, indices = [0,1,5], location = 'centroids')
         #print "2 quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values[0], [10,10,10])
-        assert allclose(quantity.vertex_values[5], [50,50,50])
+        assert numpy.allclose(quantity.vertex_values[0], [10,10,10])
+        assert numpy.allclose(quantity.vertex_values[5], [50,50,50])
         #quantity.interpolate()
         #print "quantity.centroid_values",quantity.centroid_values
-        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
+        assert numpy.allclose(quantity.centroid_values, [10,100,3,4,5,50])
 
 
@@ -2121 +2112 @@
         quantity.set_values(values, indices = [0,1,5])
         #print "quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values[0], [1,50,10])
-        assert allclose(quantity.vertex_values[5], [6,6,6])
-        assert allclose(quantity.vertex_values[1], [100,10,50])
+        assert numpy.allclose(quantity.vertex_values[0], [1,50,10])
+        assert numpy.allclose(quantity.vertex_values[5], [6,6,6])
+        assert numpy.allclose(quantity.vertex_values[1], [100,10,50])
 
         quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
@@ -2130 +2121 @@
         quantity.set_values(values, indices = [3,3,5])
         quantity.interpolate()
-        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
+        assert numpy.allclose(quantity.centroid_values, [1,2,3,400,5,999])
 
         values = [[1,1,1],[2,2,2],[3,3,3],
@@ -2142 +2133 @@
         quantity.set_values(values)
         #print "1 quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
-                                                [ 1.,  0.,  5.],
-                                                [ 5.,  6.,  1.],
-                                                [ 2.,  1.,  6.],
-                                                [ 6.,  7.,  2.],
-                                                [ 3.,  2.,  7.]])
+        assert numpy.allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
+                                                      [ 1.,  0.,  5.],
+                                                      [ 5.,  6.,  1.],
+                                                      [ 2.,  1.,  6.],
+                                                      [ 6.,  7.,  2.],
+                                                      [ 3.,  2.,  7.]])
 
     def test_setting_unique_vertex_values(self):
@@ -2155 +2146 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2171 +2161 @@
                             indices = indices)
         #print "quantity.centroid_values",quantity.centroid_values
-        assert allclose(quantity.vertex_values[0], [0,7,0])
-        assert allclose(quantity.vertex_values[1], [7,1,7])
-        assert allclose(quantity.vertex_values[2], [7,2,7])
+        assert numpy.allclose(quantity.vertex_values[0], [0,7,0])
+        assert numpy.allclose(quantity.vertex_values[1], [7,1,7])
+        assert numpy.allclose(quantity.vertex_values[2], [7,2,7])
 
 
@@ -2182 +2172 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2205 +2194 @@
 
         answer = [0.5,2,4,5,0,1,3,4.5]
-        assert allclose(answer,
-                        quantity.get_values(location = 'unique vertices'))
+        assert numpy.allclose(answer,
+                              quantity.get_values(location = 'unique vertices'))
 
         indices = [0,5,3]
         answer = [0.5,1,5]
-        assert allclose(answer,
-                        quantity.get_values(indices=indices, \
-                                            location = 'unique vertices'))
+        assert numpy.allclose(answer,
+                              quantity.get_values(indices=indices, 
+                                                  location = 'unique vertices'))
         #print "quantity.centroid_values",quantity.centroid_values
         #print "quantity.get_values(location = 'centroids') ",\
@@ -2241 +2230 @@
         quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
         
-        assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
-        assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
-
-
-        assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
-                                                                   [4,2,6],
-                                                                   [6,2,4],
-                                                                   [8,4,6]])
-        
-        assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
-                                                                                  [8,4,6]])
-
-
-        assert allclose(quantity.get_values(location='edges'), [[1,3,2],
-                                                                [4,5,3],
-                                                                [3,5,4],
-                                                                [5,7,6]])
-        assert allclose(quantity.get_values(location='edges', indices=[1,3]),
-                        [[4,5,3],
-                         [5,7,6]])        
+        assert numpy.allclose(quantity.get_values(location='centroids'), [2,4,4,6])
+        assert numpy.allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
+
+
+        assert numpy.allclose(quantity.get_values(location='vertices'), [[4,0,2],
+                                                                         [4,2,6],
+                                                                         [6,2,4],
+                                                                         [8,4,6]])
+        
+        assert numpy.allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
+                                                                                        [8,4,6]])
+
+
+        assert numpy.allclose(quantity.get_values(location='edges'), [[1,3,2],
+                                                                      [4,5,3],
+                                                                      [3,5,4],
+                                                                      [5,7,6]])
+        assert numpy.allclose(quantity.get_values(location='edges', indices=[1,3]),
+                              [[4,5,3],
+                               [5,7,6]])        
 
         # Check averaging over vertices
@@ -2280 +2269 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2298 +2286 @@
         
         #print quantity.get_values(points=interpolation_points)
-        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
+        assert numpy.allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
 
 
@@ -2311 +2299 @@
         #print answer
         #print quantity.get_values(interpolation_points=interpolation_points)
-        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points,
-                                                    verbose=False))        
+        assert numpy.allclose(answer, quantity.get_values(interpolation_points=interpolation_points,
+                                                          verbose=False))        
                         
 
@@ -2318 +2306 @@
         #indices = [0,5,3]
         #answer = [0.5,1,5]
-        #assert allclose(answer,
-        #                quantity.get_values(indices=indices, \
-        #                                    location = 'unique vertices'))
+        #assert numpy.allclose(answer,
+        #                      quantity.get_values(indices=indices,
+        #                                          location = 'unique vertices'))
 
 
@@ -2345 +2333 @@
         x, y = 2.0/3, 8.0/3
         v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 6)        
+        assert numpy.allclose(v, 6)        
 
         # Then another to test that algorithm won't blindly
@@ -2351 +2339 @@
         x, y = 4.0/3, 4.0/3
         v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 4)        
+        assert numpy.allclose(v, 4)        
 
 
@@ -2380 +2368 @@
         x, y = 2.0/3, 8.0/3
         v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]])
-        assert allclose(v, 6)
+        assert numpy.allclose(v, 6)
         
 
@@ -2387 +2375 @@
         x, y = 4.0/3, 4.0/3
         v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]])
-        assert allclose(v, 4)        
+        assert numpy.allclose(v, 4)        
         
         # Try two points
@@ -2393 +2381 @@
                [4.0/3 + xllcorner, 4.0/3 + yllcorner]]         
         v = quantity.get_values(interpolation_points=pts)
-        assert allclose(v, [6, 4])               
+        assert numpy.allclose(v, [6, 4])               
         
         # Test it using the geospatial data format with absolute input points and default georef
         pts = Geospatial_data(data_points=pts)
         v = quantity.get_values(interpolation_points=pts)
-        assert allclose(v, [6, 4])                                
+        assert numpy.allclose(v, [6, 4])                                
         
         
@@ -2405 +2393 @@
                               geo_reference=Geo_reference(zone,xllcorner,yllcorner))
         v = quantity.get_values(interpolation_points=pts)
-        assert allclose(v, [6, 4])                        
+        assert numpy.allclose(v, [6, 4])                        
         
         
@@ -2416 +2404 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2438 +2425 @@
         #print "quantity.get_values(location = 'centroids') ",\
         #      quantity.get_values(location = 'centroids')
-        assert allclose(quantity.centroid_values,
-                        quantity.get_values(location = 'centroids'))
+        assert numpy.allclose(quantity.centroid_values,
+                              quantity.get_values(location = 'centroids'))
 
 
@@ -2445 +2432 @@
         quantity.set_values(value, indices = indices)
         #print "1 quantity.vertex_values",quantity.vertex_values
-        assert allclose(quantity.vertex_values, quantity.get_values())
-
-        assert allclose(quantity.edge_values,
-                        quantity.get_values(location = 'edges'))
+        assert numpy.allclose(quantity.vertex_values, quantity.get_values())
+
+        assert numpy.allclose(quantity.edge_values,
+                              quantity.get_values(location = 'edges'))
 
         # get a subset of elements
         subset = quantity.get_values(location='centroids', indices=[0,5])
         answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
-        assert allclose(subset, answer)
+        assert numpy.allclose(subset, answer)
 
 
@@ -2460 +2447 @@
         #print "subset",subset
         #print "answer",answer
-        assert allclose(subset, answer)
+        assert numpy.allclose(subset, answer)
 
         subset = quantity.get_values( indices=[1,5])
@@ -2466 +2453 @@
         #print "subset",subset
         #print "answer",answer
-        assert allclose(subset, answer)
+        assert numpy.allclose(subset, answer)
 
     def test_smooth_vertex_values(self):
@@ -2474 +2461 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2501 +2487 @@
         
         #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5]
-        #assert allclose(answer,
-        #                quantity.get_values(location = 'unique vertices'))
+        #assert numpy.allclose(answer,
+        #                      quantity.get_values(location = 'unique vertices'))
 
         quantity.smooth_vertex_values()
@@ -2512 +2498 @@
                                 [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]]
         
-        assert allclose(answer_vertex_values,
-                        quantity.vertex_values)
+        assert numpy.allclose(answer_vertex_values,
+                              quantity.vertex_values)
         #print "quantity.centroid_values",quantity.centroid_values
         #print "quantity.get_values(location = 'centroids') ",\
@@ -2522 +2508 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Quantity, 'test')    
-    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
-
-    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
-    #print "restricted test"
-    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
+#    suite = unittest.makeSuite(Test_Quantity, 'test')    
+    suite = unittest.makeSuite(Test_Quantity, 'test_get_extrema_1')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py

@@ -6 +6 @@
 from domain import *
 from region import *
-#from anuga.config import epsilon
-from Numeric import allclose, average #, array, ones, Float
+import numpy
+
 """
 This is what the mesh in these tests look like;
@@ -44 +44 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -64 +63 @@
         domain.set_region([a, b])
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),\
-                        [[ 0.09,  0.09,  0.09],
-                         [ 0.09,  0.09,  0.09],
-                         [ 0.07,  0.07,  0.07],
-                         [ 0.07,  0.07,  0.07],
-                         [ 1.0,  1.0,  1.0],
-                         [ 1.0,  1.0,  1.0]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),
+                              [[ 0.09,  0.09,  0.09],
+                               [ 0.09,  0.09,  0.09],
+                               [ 0.07,  0.07,  0.07],
+                               [ 0.07,  0.07,  0.07],
+                               [ 1.0,  1.0,  1.0],
+                               [ 1.0,  1.0,  1.0]])
 
         #c = Add_Value_To_region('all', 'friction', 10.0)
         domain.set_region(Add_value_to_region('all', 'friction', 10.0))
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),
-                        [[ 10.09, 10.09, 10.09],
-                         [ 10.09, 10.09, 10.09],
-                         [ 10.07, 10.07, 10.07],
-                         [ 10.07, 10.07, 10.07],
-                         [ 11.0,  11.0,  11.0],
-                         [ 11.0,  11.0,  11.0]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),
+                              [[ 10.09, 10.09, 10.09],
+                               [ 10.09, 10.09, 10.09],
+                               [ 10.07, 10.07, 10.07],
+                               [ 10.07, 10.07, 10.07],
+                               [ 11.0,  11.0,  11.0],
+                               [ 11.0,  11.0,  11.0]])
 
         # trying a function
         domain.set_region(Set_region('top', 'friction', add_x_y))
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),
-                        [[ 10.09, 10.09, 10.09],
-                         [ 10.09, 10.09, 10.09],
-                         [ 10.07, 10.07, 10.07],
-                         [ 10.07, 10.07, 10.07],
-                         [ 5./3,  2.0,  2./3],
-                         [ 1.0,  2./3,  2.0]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),
+                              [[ 10.09, 10.09, 10.09],
+                               [ 10.09, 10.09, 10.09],
+                               [ 10.07, 10.07, 10.07],
+                               [ 10.07, 10.07, 10.07],
+                               [ 5./3,  2.0,  2./3],
+                               [ 1.0,  2./3,  2.0]])
 
         domain.set_quantity('elevation', 10.0)
@@ -98 +97 @@
         domain.set_region(Add_value_to_region('top', 'stage', 1.0,initial_quantity='elevation'))
         #print domain.quantities['stage'].get_values()
-        assert allclose(domain.quantities['stage'].get_values(),
-                        [[ 10., 10., 10.],
-                         [ 10., 10., 10.],
-                         [ 10., 10., 10.],
-                         [ 10., 10., 10.],
-                         [ 11.0,  11.0,  11.0],
-                         [ 11.0,  11.0,  11.0]])
+        assert numpy.allclose(domain.quantities['stage'].get_values(),
+                              [[ 10., 10., 10.],
+                               [ 10., 10., 10.],
+                               [ 10., 10., 10.],
+                               [ 10., 10., 10.],
+                               [ 11.0,  11.0,  11.0],
+                               [ 11.0,  11.0,  11.0]])
 
         
@@ -115 +114 @@
         domain.set_region(Add_quantities('top', 'elevation','stage'))
         #print domain.quantities['stage'].get_values()
-        assert allclose(domain.quantities['elevation'].get_values(),
-                        [[ 10., 10., 10.],
-                         [ 10., 10., 10.],
-                         [ 10., 10., 10.],
-                         [ 10., 10., 10.],
-                         [ 33.,  33.0,  33.],
-                         [ 33.0,  33.,  33.]])
+        assert numpy.allclose(domain.quantities['elevation'].get_values(),
+                              [[ 10., 10., 10.],
+                               [ 10., 10., 10.],
+                               [ 10., 10., 10.],
+                               [ 10., 10., 10.],
+                               [ 33.,  33.0,  33.],
+                               [ 33.0,  33.,  33.]])
         
     def test_unique_vertices(self):
@@ -129 +128 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -147 +145 @@
         domain.set_region(a)
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),\
-                        [[ 0.09,  0.09,  0.09],
-                         [ 0.09,  0.09,  0.09],
-                         [ 0.09,  0.07,  0.09],
-                         [ 0.07,  0.09,  0.07],
-                         [ 0.07,  0.07,  0.07],
-                         [ 0.07,  0.07,  0.07]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),\
+                              [[ 0.09,  0.09,  0.09],
+                               [ 0.09,  0.09,  0.09],
+                               [ 0.09,  0.07,  0.09],
+                               [ 0.07,  0.09,  0.07],
+                               [ 0.07,  0.07,  0.07],
+                               [ 0.07,  0.07,  0.07]])
 
 
@@ -162 +160 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -180 +177 @@
 
         #print domain.quantities['friction'].get_values()
-        assert allclose(domain.quantities['friction'].get_values(),\
-                        [[ 1.07,  1.07,  1.07],
-                         [ 1.07,  1.07,  1.07],
-                         [ 1.07,  0.07,  1.07],
-                         [ 0.07,  1.07,  0.07],
-                         [ 0.07,  0.07,  0.07],
-                         [ 0.07,  0.07,  0.07]])
+        assert numpy.allclose(domain.quantities['friction'].get_values(),\
+                              [[ 1.07,  1.07,  1.07],
+                               [ 1.07,  1.07,  1.07],
+                               [ 1.07,  0.07,  1.07],
+                               [ 0.07,  1.07,  0.07],
+                               [ 0.07,  0.07,  0.07],
+                               [ 0.07,  0.07,  0.07]])
                          
     def test_unique_vertices_average_loc_vert(self):
@@ -194 +191 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -218 +214 @@
         #print domain.quantities['friction'].get_values()
         frict_points = domain.quantities['friction'].get_values()
-        assert allclose(frict_points[0],\
-                        [ calc_frict, calc_frict, calc_frict])
-        assert allclose(frict_points[1],\
-                        [ calc_frict, calc_frict, calc_frict])
+        assert numpy.allclose(frict_points[0],
+                              [ calc_frict, calc_frict, calc_frict])
+        assert numpy.allclose(frict_points[1],
+                              [ calc_frict, calc_frict, calc_frict])
   
     def test_unique_vertices_average_loc_unique_vert(self):
@@ -229 +225 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -252 +247 @@
         #print domain.quantities['friction'].get_values()
         frict_points = domain.quantities['friction'].get_values()
-        assert allclose(frict_points[0],\
-                        [ calc_frict, calc_frict, calc_frict])
-        assert allclose(frict_points[1],\
-                        [ calc_frict, calc_frict, calc_frict])
-        assert allclose(frict_points[2],\
-                        [ calc_frict, 1.0 + 2.0/3.0, calc_frict])
-        assert allclose(frict_points[3],\
-                        [ 2.0/3.0,calc_frict, 1.0 + 2.0/3.0])
+        assert numpy.allclose(frict_points[0],
+                              [ calc_frict, calc_frict, calc_frict])
+        assert numpy.allclose(frict_points[1],
+                              [ calc_frict, calc_frict, calc_frict])
+        assert numpy.allclose(frict_points[2],
+                              [ calc_frict, 1.0 + 2.0/3.0, calc_frict])
+        assert numpy.allclose(frict_points[3],
+                              [ 2.0/3.0,calc_frict, 1.0 + 2.0/3.0])
                                                 
                          
@@ -265 +260 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Region,'test')
+##    suite = unittest.makeSuite(Test_Region,'test_unique_vertices_average_loc_unique_vert')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

@@ -1 +1 @@
 #!/usr/bin/env python
 
-
 import unittest
-from Numeric import zeros, array, allclose, Float
+import numpy
+
 from math import sqrt, pi
 import tempfile, os
@@ -91 +91 @@
 
             #Exact linear intpolation
-            assert allclose(q[0], 2*t)
+            assert numpy.allclose(q[0], 2*t)
             if i%6 == 0:
-                assert allclose(q[1], t**2)
-                assert allclose(q[2], sin(t*pi/600))
+                assert numpy.allclose(q[1], t**2)
+                assert numpy.allclose(q[2], sin(t*pi/600))
 
         #Check non-exact
@@ -100 +100 @@
         t = 90 #Halfway between 60 and 120
         q = F(t)
-        assert allclose( (120**2 + 60**2)/2, q[1] )
-        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
+        assert numpy.allclose( (120**2 + 60**2)/2, q[1] )
+        assert numpy.allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
 
 
         t = 100 #Two thirds of the way between between 60 and 120
         q = F(t)
-        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
-        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
+        assert numpy.allclose( 2*120**2/3 + 60**2/3, q[1] )
+        assert numpy.allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
 
         os.remove(filename + '.txt')
@@ -125 +125 @@
         from shallow_water import Domain, Dirichlet_boundary
         from mesh_factory import rectangular
-        from Numeric import take, concatenate, reshape
 
         #Create basic mesh and shallow water domain
@@ -148 +147 @@
 
         # Boundary conditions
+        print 'test_util: 0'
         B0 = Dirichlet_boundary([0,0,0])
         B6 = Dirichlet_boundary([0.6,0,0])
+        print 'test_util: 1'
         domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
         domain1.check_integrity()
+        print 'test_util: 2'
 
         finaltime = 8
         #Evolution
         t0 = -1
+        print 'test_util: 3'
+# crash at following line
         for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
+            print 'test_util: 4'
             #print 'Timesteps: %.16f, %.16f' %(t0, t)
             #if t == t0:
@@ -164 +169 @@
              
             #domain1.write_time()
+        print 'test_util: 5'
 
 
@@ -182 +188 @@
 
         last_time_index = len(time)-1 #Last last_time_index
-        d_stage = reshape(take(stage[last_time_index, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(numpy.take(stage[last_time_index, :], [0,5,10,15], axis=0), (4,1))
+        d_uh = numpy.reshape(numpy.take(xmomentum[last_time_index, :], [0,5,10,15], axis=0), (4,1))
+        d_vh = numpy.reshape(numpy.take(ymomentum[last_time_index, :], [0,5,10,15], axis=0), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -194 +200 @@
 
         #And the midpoints are found now
-        Dx = take(reshape(x, (16,1)), [0,5,10,15])
-        Dy = take(reshape(y, (16,1)), [0,5,10,15])
-
-        diag = concatenate( (Dx, Dy), axis=1)
+        Dx = numpy.take(numpy.reshape(x, (16,1)), [0,5,10,15], axis=0)
+        Dy = numpy.take(numpy.reshape(y, (16,1)), [0,5,10,15], axis=0)
+
+        diag = numpy.concatenate( (Dx, Dy), axis=1)
         d_midpoints = (diag[1:] + diag[:-1])/2
 
@@ -209 +215 @@
         assert not T[-1] == T[-2], msg
         t = time[last_time_index]
-        q = f(t, point_id=0); assert allclose(r0, q)
-        q = f(t, point_id=1); assert allclose(r1, q)
-        q = f(t, point_id=2); assert allclose(r2, q)
+        q = f(t, point_id=0); assert numpy.allclose(r0, q)
+        q = f(t, point_id=1); assert numpy.allclose(r1, q)
+        q = f(t, point_id=2); assert numpy.allclose(r2, q)
 
 
@@ -218 +224 @@
 
         timestep = 0 #First timestep
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(numpy.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
+        d_uh = numpy.reshape(numpy.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
+        d_vh = numpy.reshape(numpy.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -231 +237 @@
         #Let us see if the file function can find the correct
         #values
-        q = f(0, point_id=0); assert allclose(r0, q)
-        q = f(0, point_id=1); assert allclose(r1, q)
-        q = f(0, point_id=2); assert allclose(r2, q)
+        q = f(0, point_id=0); assert numpy.allclose(r0, q)
+        q = f(0, point_id=1); assert numpy.allclose(r1, q)
+        q = f(0, point_id=2); assert numpy.allclose(r2, q)
 
 
@@ -240 +246 @@
 
         timestep = 33
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(numpy.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
+        d_uh = numpy.reshape(numpy.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
+        d_vh = numpy.reshape(numpy.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -251 +257 @@
         r2 = (D[2] + D[3])/2
 
-        q = f(timestep/10., point_id=0); assert allclose(r0, q)
-        q = f(timestep/10., point_id=1); assert allclose(r1, q)
-        q = f(timestep/10., point_id=2); assert allclose(r2, q)
+        q = f(timestep/10., point_id=0); assert numpy.allclose(r0, q)
+        q = f(timestep/10., point_id=1); assert numpy.allclose(r1, q)
+        q = f(timestep/10., point_id=2); assert numpy.allclose(r2, q)
 
 
@@ -261 +267 @@
 
         timestep = 15
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(numpy.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
+        d_uh = numpy.reshape(numpy.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
+        d_vh = numpy.reshape(numpy.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -274 +280 @@
         #
         timestep = 16
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(numpy.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
+        d_uh = numpy.reshape(numpy.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
+        d_vh = numpy.reshape(numpy.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -290 +296 @@
         r2 = (r2_0 + r2_1)/2
 
-        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
-        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
-        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
+        q = f((timestep - 0.5)/10., point_id=0); assert numpy.allclose(r0, q)
+        q = f((timestep - 0.5)/10., point_id=1); assert numpy.allclose(r1, q)
+        q = f((timestep - 0.5)/10., point_id=2); assert numpy.allclose(r2, q)
 
         ##################
@@ -304 +310 @@
 
         #And the file function gives
-        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
-        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
-        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
+        q = f((timestep - 1.0/3)/10., point_id=0); assert numpy.allclose(r0, q)
+        q = f((timestep - 1.0/3)/10., point_id=1); assert numpy.allclose(r1, q)
+        q = f((timestep - 1.0/3)/10., point_id=2); assert numpy.allclose(r2, q)
 
         fid.close()
@@ -326 +332 @@
         from shallow_water import Domain, Dirichlet_boundary
         from mesh_factory import rectangular
-        from Numeric import take, concatenate, reshape
 
 
@@ -386 +391 @@
 
         last_time_index = len(time)-1 #Last last_time_index     
-        d_stage = reshape(take(stage[last_time_index, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(take(stage[last_time_index, :], [0,5,10,15]), (4,1))
+        d_uh = numpy.reshape(take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))
+        d_vh = numpy.reshape(take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -398 +403 @@
 
         #And the midpoints are found now
-        Dx = take(reshape(x, (16,1)), [0,5,10,15])
-        Dy = take(reshape(y, (16,1)), [0,5,10,15])
-
-        diag = concatenate( (Dx, Dy), axis=1)
+        Dx = take(numpy.reshape(x, (16,1)), [0,5,10,15])
+        Dy = take(numpy.reshape(y, (16,1)), [0,5,10,15])
+
+        diag = numpy.concatenate( (Dx, Dy), axis=1)
         d_midpoints = (diag[1:] + diag[:-1])/2
 
@@ -415 +420 @@
 
         t = time[last_time_index]                         
-        q = f(t, point_id=0); assert allclose(r0, q)
-        q = f(t, point_id=1); assert allclose(r1, q)
-        q = f(t, point_id=2); assert allclose(r2, q)
+        q = f(t, point_id=0); assert numpy.allclose(r0, q)
+        q = f(t, point_id=1); assert numpy.allclose(r1, q)
+        q = f(t, point_id=2); assert numpy.allclose(r2, q)
 
 
@@ -424 +429 @@
 
         timestep = 0 #First timestep
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = numpy.reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = numpy.reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -437 +442 @@
         #Let us see if the file function can find the correct
         #values
-        q = f(0, point_id=0); assert allclose(r0, q)
-        q = f(0, point_id=1); assert allclose(r1, q)
-        q = f(0, point_id=2); assert allclose(r2, q)
+        q = f(0, point_id=0); assert numpy.allclose(r0, q)
+        q = f(0, point_id=1); assert numpy.allclose(r1, q)
+        q = f(0, point_id=2); assert numpy.allclose(r2, q)
 
 
@@ -446 +451 @@
 
         timestep = 33
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = numpy.reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = numpy.reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -457 +462 @@
         r2 = (D[2] + D[3])/2
 
-        q = f(timestep/10., point_id=0); assert allclose(r0, q)
-        q = f(timestep/10., point_id=1); assert allclose(r1, q)
-        q = f(timestep/10., point_id=2); assert allclose(r2, q)
+        q = f(timestep/10., point_id=0); assert numpy.allclose(r0, q)
+        q = f(timestep/10., point_id=1); assert numpy.allclose(r1, q)
+        q = f(timestep/10., point_id=2); assert numpy.allclose(r2, q)
 
 
@@ -467 +472 @@
 
         timestep = 15
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = numpy.reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = numpy.reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -480 +485 @@
         #
         timestep = 16
-        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
-        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
-        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
-        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+        d_stage = numpy.reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = numpy.reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = numpy.reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = numpy.concatenate( (d_stage, d_uh, d_vh), axis=1)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -496 +501 @@
         r2 = (r2_0 + r2_1)/2
 
-        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
-        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
-        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
+        q = f((timestep - 0.5)/10., point_id=0); assert numpy.allclose(r0, q)
+        q = f((timestep - 0.5)/10., point_id=1); assert numpy.allclose(r1, q)
+        q = f((timestep - 0.5)/10., point_id=2); assert numpy.allclose(r2, q)
 
         ##################
@@ -510 +515 @@
 
         #And the file function gives
-        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
-        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
-        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
+        q = f((timestep - 1.0/3)/10., point_id=0); assert numpy.allclose(r0, q)
+        q = f((timestep - 1.0/3)/10., point_id=1); assert numpy.allclose(r1, q)
+        q = f((timestep - 1.0/3)/10., point_id=2); assert numpy.allclose(r2, q)
 
         fid.close()
@@ -542 +547 @@
         import os, time
         from anuga.config import time_format
-        from Numeric import sin, pi, exp
         from mesh_factory import rectangular
         from shallow_water import Domain
@@ -584 +588 @@
             domain.set_quantity('xmomentum', f2)
 
-            f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
+            f3 = lambda x,y: x**2 + y**2 * numpy.sin(t*numpy.pi/600)
             domain.set_quantity('ymomentum', f3)
 
@@ -604 +608 @@
 
         #Check that FF updates fixes domain starttime
-        assert allclose(domain.starttime, start)
+        assert numpy.allclose(domain.starttime, start)
 
         #Check that domain.starttime isn't updated if later
@@ -611 +615 @@
                           quantities = domain.conserved_quantities,
                           interpolation_points = interpolation_points)
-        assert allclose(domain.starttime, start+1)
+        assert numpy.allclose(domain.starttime, start+1)
         domain.starttime = start
 
@@ -645 +649 @@
                      self.failUnless( q == actual, 'Fail!')
                 else:
-                    assert allclose(q, actual)
+                    assert numpy.allclose(q, actual)
 
 
@@ -655 +659 @@
             t = 90 #Halfway between 60 and 120
             q = F(t, point_id=id)
-            assert allclose( (q120+q60)/2, q )
+            assert numpy.allclose( (q120+q60)/2, q )
 
             t = 100 #Two thirds of the way between between 60 and 120
             q = F(t, point_id=id)
-            assert allclose(q60/3 + 2*q120/3, q)
+            assert numpy.allclose(q60/3 + 2*q120/3, q)
 
 
@@ -670 +674 @@
                           quantities = domain.conserved_quantities,
                           interpolation_points = interpolation_points)
-        assert allclose(domain.starttime, start+delta)
+        assert numpy.allclose(domain.starttime, start+delta)
 
 
@@ -689 +693 @@
 
                 q = F(t-delta, point_id=id)
-                assert allclose(q, (k*q1 + (6-k)*q0)/6)
+                assert numpy.allclose(q, (k*q1 + (6-k)*q0)/6)
 
 
@@ -703 +707 @@
         import os, time
         from anuga.config import time_format
-        from Numeric import sin, pi, exp
         from mesh_factory import rectangular
         from shallow_water import Domain
@@ -751 +754 @@
             domain.set_quantity('xmomentum', f2)
 
-            f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
+            f3 = lambda x,y: x**2 + y**2 * numpy.sin(t*numpy.pi/600)
             domain.set_quantity('ymomentum', f3)
 
@@ -774 +777 @@
 
         #Check that FF updates fixes domain starttime
-        assert allclose(domain.starttime, start)
+        assert numpy.allclose(domain.starttime, start)
 
         #Check that domain.starttime isn't updated if later
@@ -781 +784 @@
                           quantities = domain.conserved_quantities,
                           interpolation_points = interpolation_points)
-        assert allclose(domain.starttime, start+1)
+        assert numpy.allclose(domain.starttime, start+1)
         domain.starttime = start
 
@@ -817 +820 @@
                      self.failUnless( q == actual, 'Fail!')
                 else:
-                    assert allclose(q, actual)
+                    assert numpy.allclose(q, actual)
 
         # now lets check points inside the mesh
@@ -863 +866 @@
                      self.failUnless( q == actual, 'Fail!')
                 else:
-                    assert allclose(q, actual)
+                    assert numpy.allclose(q, actual)
 
 
@@ -873 +876 @@
             t = 90 #Halfway between 60 and 120
             q = F(t, point_id=id)
-            assert allclose( (q120+q60)/2, q )
+            assert numpy.allclose( (q120+q60)/2, q )
 
             t = 100 #Two thirds of the way between between 60 and 120
             q = F(t, point_id=id)
-            assert allclose(q60/3 + 2*q120/3, q)
+            assert numpy.allclose(q60/3 + 2*q120/3, q)
 
 
@@ -888 +891 @@
                           quantities = domain.conserved_quantities,
                           interpolation_points = interpolation_points)
-        assert allclose(domain.starttime, start+delta)
+        assert numpy.allclose(domain.starttime, start+delta)
 
 
@@ -907 +910 @@
 
                 q = F(t-delta, point_id=id)
-                assert allclose(q, (k*q1 + (6-k)*q0)/6)
+                assert numpy.allclose(q, (k*q1 + (6-k)*q0)/6)
 
 
@@ -955 +958 @@
                           domain,
                           quantities = ['Attribute0', 'Attribute1', 'Attribute2'])  
-        assert allclose(domain.starttime, start)
+        assert numpy.allclose(domain.starttime, start)
 
         # Check that domain.starttime is updated if too early
@@ -962 +965 @@
                           domain,
                           quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
-        assert allclose(domain.starttime, start)
+        assert numpy.allclose(domain.starttime, start)
 
         # Check that domain.starttime isn't updated if later
@@ -969 +972 @@
                           domain,
                           quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
-        assert allclose(domain.starttime, start+1)
+        assert numpy.allclose(domain.starttime, start+1)
 
         domain.starttime = start
@@ -987 +990 @@
 
             #Exact linear intpolation
-            assert allclose(q[0], 2*t)
+            assert numpy.allclose(q[0], 2*t)
             if i%6 == 0:
-                assert allclose(q[1], t**2)
-                assert allclose(q[2], sin(t*pi/600))
+                assert numpy.allclose(q[1], t**2)
+                assert numpy.allclose(q[2], sin(t*pi/600))
 
         #Check non-exact
@@ -996 +999 @@
         t = 90 #Halfway between 60 and 120
         q = F(t)
-        assert allclose( (120**2 + 60**2)/2, q[1] )
-        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
+        assert numpy.allclose( (120**2 + 60**2)/2, q[1] )
+        assert numpy.allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
 
 
         t = 100 #Two thirds of the way between between 60 and 120
         q = F(t)
-        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
-        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
+        assert numpy.allclose( 2*120**2/3 + 60**2/3, q[1] )
+        assert numpy.allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
 
         os.remove(filename + '.tms')
@@ -1052 +1055 @@
         F = file_function(filename + '.tms', domain,
                           quantities = ['Attribute0', 'Attribute1', 'Attribute2'])        
-        assert allclose(domain.starttime, start+delta)
+        assert numpy.allclose(domain.starttime, start+delta)
 
 
@@ -1063 +1066 @@
 
             #Exact linear intpolation
-            assert allclose(q[0], 2*t)
+            assert numpy.allclose(q[0], 2*t)
             if i%6 == 0:
-                assert allclose(q[1], t**2)
-                assert allclose(q[2], sin(t*pi/600))
+                assert numpy.allclose(q[1], t**2)
+                assert numpy.allclose(q[2], sin(t*pi/600))
 
         #Check non-exact
@@ -1072 +1075 @@
         t = 90 #Halfway between 60 and 120
         q = F(t-delta)
-        assert allclose( (120**2 + 60**2)/2, q[1] )
-        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
+        assert numpy.allclose( (120**2 + 60**2)/2, q[1] )
+        assert numpy.allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
 
 
         t = 100 #Two thirds of the way between between 60 and 120
         q = F(t-delta)
-        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
-        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
+        assert numpy.allclose( 2*120**2/3 + 60**2/3, q[1] )
+        assert numpy.allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
 
 
@@ -1091 +1094 @@
         #FIXME: Division is not expected to work for integers.
         #This must be caught.
-        foo = array([[1,2,3],
-                     [4,5,6]], Float)
-
-        bar = array([[-1,0,5],
-                     [6,1,1]], Float)                  
+        foo = numpy.array([[1,2,3], [4,5,6]], numpy.float)
+
+        bar = numpy.array([[-1,0,5], [6,1,1]], numpy.float)                  
 
         D = {'X': foo, 'Y': bar}
 
         Z = apply_expression_to_dictionary('X+Y', D)        
-        assert allclose(Z, foo+bar)
+        assert numpy.allclose(Z, foo+bar)
 
         Z = apply_expression_to_dictionary('X*Y', D)        
-        assert allclose(Z, foo*bar)        
+        assert numpy.allclose(Z, foo*bar)        
 
         Z = apply_expression_to_dictionary('4*X+Y', D)        
-        assert allclose(Z, 4*foo+bar)        
+        assert numpy.allclose(Z, 4*foo+bar)        
 
         # test zero division is OK
         Z = apply_expression_to_dictionary('X/Y', D)
-        assert allclose(1/Z, 1/(foo/bar)) # can't compare inf to inf
+        assert numpy.allclose(1/Z, 1/(foo/bar)) # can't compare inf to inf
 
         # make an error for zero on zero
@@ -1268 +1269 @@
         tris = [[0,1,2]]
         new_verts, new_tris = remove_lone_verts(verts, tris)
-        assert new_verts == verts
-        assert new_tris == tris
+        assert numpy.alltrue(new_verts == verts)
+        assert numpy.alltrue(new_tris == tris)
      
 
@@ -1284 +1285 @@
         new_verts, new_tris = remove_lone_verts(verts, tris)
         #print "new_verts", new_verts
-        assert new_verts == [[0,0],[1,0],[0,1]]
-        assert new_tris == [[0,1,2]]
+        assert numpy.alltrue(new_verts == [[0,0],[1,0],[0,1]])
+        assert numpy.alltrue(new_tris == [[0,1,2]])
      
     def test_remove_lone_verts_c(self):
@@ -1291 +1292 @@
         tris = [[0,1,3]]
         new_verts, new_tris = remove_lone_verts(verts, tris)
-        #print "new_verts", new_verts
-        assert new_verts == [[0,0],[1,0],[0,1]]
-        assert new_tris == [[0,1,2]]
+        print "new_verts", new_verts
+        assert numpy.alltrue(new_verts == [[0,0],[1,0],[0,1]])
+        assert numpy.alltrue(new_tris == [[0,1,2]])
         
     def test_remove_lone_verts_b(self):
@@ -1299 +1300 @@
         tris = [[0,1,2]]
         new_verts, new_tris = remove_lone_verts(verts, tris)
-        assert new_verts == verts[0:3]
-        assert new_tris == tris
+        assert numpy.alltrue(new_verts == verts[0:3])
+        assert numpy.alltrue(new_tris == tris)
      
 
@@ -1307 +1308 @@
         tris = [[0,1,2]]
         new_verts, new_tris = remove_lone_verts(verts, tris)
-        assert new_verts == verts[0:3]
-        assert new_tris == tris
+        assert numpy.alltrue(new_verts == verts[0:3])
+        assert numpy.alltrue(new_tris == tris)
         
     def test_get_min_max_values(self):
@@ -1496 +1497 @@
             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
             #print 'assert line',line[i],'point1',point1_answers_array[i]
-            assert allclose(line[i], point1_answers_array[i])
+            assert numpy.allclose(line[i], point1_answers_array[i])
 
         point2_answers_array = [[0.0,1.0,1.5,-0.5,3.0,4.0], [2.0,10.0,10.5,-0.5,3.0,4.0]]
@@ -1509 +1510 @@
             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
             #print 'assert line',line[i],'point1',point1_answers_array[i]
-            assert allclose(line[i], point2_answers_array[i])
+            assert numpy.allclose(line[i], point2_answers_array[i])
                          
         # clean up
@@ -1615 +1616 @@
             line.append([float(row[0]),float(row[1]),float(row[2])])
             #print 'line',line[i],'point1',point1_answers_array[i]
-            assert allclose(line[i], point1_answers_array[i])
+            assert numpy.allclose(line[i], point1_answers_array[i])
 
         point2_answers_array = [[0.0,1.0,-0.5], [2.0,10.0,-0.5]]
@@ -1628 +1629 @@
             line.append([float(row[0]),float(row[1]),float(row[2])])
 #            print 'line',line[i],'point1',point1_answers_array[i]
-            assert allclose(line[i], point2_answers_array[i])
+            assert numpy.allclose(line[i], point2_answers_array[i])
                          
         # clean up
@@ -1734 +1735 @@
             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
             #print 'assert line',line[i],'point1',point1_answers_array[i]
-            assert allclose(line[i], point1_answers_array[i])
+            assert numpy.allclose(line[i], point1_answers_array[i])
 
         point2_answers_array = [[5.0,1.0,1.5,-0.5,3.0,4.0], [7.0,10.0,10.5,-0.5,3.0,4.0]]
@@ -1747 +1748 @@
             line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4]),float(row[5])])
             #print 'assert line',line[i],'point1',point1_answers_array[i]
-            assert allclose(line[i], point2_answers_array[i])
+            assert numpy.allclose(line[i], point2_answers_array[i])
                          
         # clean up
@@ -1830 +1831 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Util,'test')
-#    suite = unittest.makeSuite(Test_Util,'test_sww2csv')
+#    suite = unittest.makeSuite(Test_Util,'test_spatio_temporal_file_function_basic')
 #    runner = unittest.TextTestRunner(verbosity=2)
     runner = unittest.TextTestRunner(verbosity=1)

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -15 +15 @@
 
 from anuga.utilities.numerical_tools import ensure_numeric
-from Numeric import arange, choose, zeros, Float, array, allclose, take, compress
+import numpy
     
 from anuga.geospatial_data.geospatial_data import ensure_absolute
@@ -237 +237 @@
     from anuga.config import time_format
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
+##    from numpy.oldnumeric import array, zeros, Float, alltrue, concatenate, reshape
 
     # Open NetCDF file
@@ -303 +303 @@
     # Get variables
     # if verbose: print 'Get variables'    
-    time = fid.variables['time'][:]    
+    time = numpy.array(fid.variables['time'][:]    )
 
     # Get time independent stuff
@@ -317 +317 @@
             triangles = fid.variables['volumes'][:]
 
-        x = reshape(x, (len(x),1))
-        y = reshape(y, (len(y),1))
-        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
+        x = numpy.reshape(x, (len(x),1))
+        y = numpy.reshape(y, (len(y),1))
+        vertex_coordinates = numpy.concatenate((x,y), axis=1) #m x 2 array
 
         if boundary_polygon is not None:
@@ -331 +331 @@
             for i in range(len(boundary_polygon)):
                 for j in range(len(x)):
-                    if allclose(vertex_coordinates[j],boundary_polygon[i],1e-4):
+                    if numpy.allclose(vertex_coordinates[j],boundary_polygon[i],1e-4):
                         #FIX ME:
                         #currently gauges lat and long is stored as float and
@@ -351 +351 @@
                 gauge_neighbour_id.append(-1)
             gauge_neighbour_id=ensure_numeric(gauge_neighbour_id)
-            if len(compress(gauge_neighbour_id>=0,gauge_neighbour_id))!=len(temp)-1:
+            if len(numpy.compress(gauge_neighbour_id>=0,gauge_neighbour_id))!=len(temp)-1:
                 msg='incorrect number of segments'
                 raise msg
@@ -373 +373 @@
         # move time back - relative to domain's time
         if domain_starttime > starttime:
+##            print 'type(time)=%s' % type(time)
+##            print 'type(domain_starttime)=%s' % type(domain_starttime)
+##            print 'type(starttime)=%s' % type(starttime)
+##            a = time - domain_starttime
             time = time - domain_starttime + starttime
 
@@ -398 +402 @@
         if boundary_polygon is not None:
             #removes sts points that do not lie on boundary
-            quantities[name] = take(quantities[name],gauge_id,1)
+            quantities[name] = numpy.take(quantities[name], gauge_id, 1)
             
     # Close sww, tms or sts netcdf file         
@@ -520 +524 @@
                 raise 'Illegal input to get_textual_float:', value
         else:
-            return format %float(value)
+            return format % float(value)
 
 
@@ -1036 +1040 @@
     """
 #    from math import sqrt, atan, degrees
-    from Numeric import ones, allclose, zeros, Float, ravel
+##    from numpy.oldnumeric import ones, allclose, zeros, Float, ravel
     from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
 
@@ -1077 +1081 @@
     n0 = int(n0)
     m = len(locations)
-    model_time = zeros((n0,m,p), Float) 
-    stages = zeros((n0,m,p), Float)
-    elevations = zeros((n0,m,p), Float) 
-    momenta = zeros((n0,m,p), Float)
-    xmom = zeros((n0,m,p), Float)
-    ymom = zeros((n0,m,p), Float)
-    speed = zeros((n0,m,p), Float)
-    bearings = zeros((n0,m,p), Float)
-    due_east = 90.0*ones((n0,1), Float)
-    due_west = 270.0*ones((n0,1), Float)
-    depths = zeros((n0,m,p), Float)
-    eastings = zeros((n0,m,p), Float)
+    model_time = numpy.zeros((n0,m,p), numpy.float) 
+    stages = numpy.zeros((n0,m,p), numpy.float)
+    elevations = numpy.zeros((n0,m,p), numpy.float) 
+    momenta = numpy.zeros((n0,m,p), numpy.float)
+    xmom = numpy.zeros((n0,m,p), numpy.float)
+    ymom = numpy.zeros((n0,m,p), numpy.float)
+    speed = numpy.zeros((n0,m,p), numpy.float)
+    bearings = numpy.zeros((n0,m,p), numpy.float)
+    due_east = 90.0*numpy.ones((n0,1), numpy.float)
+    due_west = 270.0*numpy.ones((n0,1), numpy.float)
+    depths = numpy.zeros((n0,m,p), numpy.float)
+    eastings = numpy.zeros((n0,m,p), numpy.float)
     min_stages = []
     max_stages = []
@@ -1100 +1104 @@
     min_speeds = []    
     max_depths = []
-    model_time_plot3d = zeros((n0,m), Float)
-    stages_plot3d = zeros((n0,m), Float)
-    eastings_plot3d = zeros((n0,m),Float)
+    model_time_plot3d = numpy.zeros((n0,m), numpy.float)
+    stages_plot3d = numpy.zeros((n0,m), numpy.float)
+    eastings_plot3d = numpy.zeros((n0,m),numpy.float)
     if time_unit is 'mins': scale = 60.0
     if time_unit is 'hours': scale = 3600.0
@@ -1202 +1206 @@
                 else:
                     #ax.plot_wireframe(model_time[:,:,j],eastings[:,:,j],stages[:,:,j])
-                    ax.plot3D(ravel(eastings[:,:,j]),ravel(model_time[:,:,j]),ravel(stages[:,:,j]))
+                    ax.plot3D(numpy.ravel(eastings[:,:,j]),numpy.ravel(model_time[:,:,j]),numpy.ravel(stages[:,:,j]))
                 ax.set_xlabel('time')
                 ax.set_ylabel('x')
@@ -1501 +1505 @@
     
     # initialise the array to easily find the index of the first loner
-    loners=arange(2*N, N, -1) # if N=3 [6,5,4]
+    loners=numpy.arange(2*N, N, -1) # if N=3 [6,5,4]
     for t in triangles:
         for vert in t:
@@ -1535 +1539 @@
         #print "loners", loners
         #print "triangles before", triangles
-        triangles = choose(triangles,loners)
+        triangles = numpy.choose(triangles,loners)
         #print "triangles after", triangles
     return verts, triangles
@@ -1550 +1554 @@
 
         
-    xc = zeros(triangles.shape[0], Float) # Space for centroid info
+    xc = numpy.zeros(triangles.shape[0], numpy.float) # Space for centroid info
     
     for k in range(triangles.shape[0]):
@@ -1832 +1836 @@
                 #add tide to stage if provided
                 if quantity == 'stage':
-                     quantity_value[quantity]=array(quantity_value[quantity])+directory_add_tide
+                     quantity_value[quantity]=numpy.array(quantity_value[quantity])+directory_add_tide
 
                 #condition to find max and mins for all the plots
@@ -1943 +1947 @@
             #get data from dict in to list
             #do maths to list by changing to array
-            t=(array(directory_quantity_value[directory][filename]['time'])+directory_start_time)/seconds_in_minutes
+            t=(numpy.array(directory_quantity_value[directory][filename]['time'])+directory_start_time)/seconds_in_minutes
 
             #finds the maximum elevation, used only as a test
@@ -2124 +2128 @@
     from csv import reader,writer
     from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN
-    from Numeric import array, resize, shape, Float, zeros, take, argsort, argmin
+##    from numpy.oldnumeric import array, resize, shape, Float, zeros, take, argsort, argmin
     import string
     from anuga.shallow_water.data_manager import get_all_swwfiles
@@ -2165 +2169 @@
         
     #convert to array for file_function
-    points_array = array(points,Float)
+    points_array = numpy.array(points, numpy.float)
         
     points_array = ensure_absolute(points_array)