Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5892)
@@ -8,5 +8,7 @@
 """
 
-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,5 +107,4 @@
 
         if verbose: print 'Initialising Domain'
-        from Numeric import zeros, Float, Int, ones
         from quantity import Quantity
 
@@ -155,10 +156,10 @@
         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,5 +169,5 @@
         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,5 +176,5 @@
         # 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,10 +231,10 @@
         # 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,6 +267,4 @@
         """
 
-        from Numeric import zeros, Float
-
         if not (vertex is None or edge is None):
             msg = 'Values for both vertex and edge was specified.'
@@ -273,5 +272,5 @@
             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,5 +768,5 @@
             # 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,4 +1085,5 @@
 
         """
+        print '.evolve: 0'
 
         from anuga.config import min_timestep, max_timestep, epsilon
@@ -1096,4 +1096,5 @@
                %self.get_boundary_tags()
         assert hasattr(self, 'boundary_objects'), msg
+        print '.evolve: 1'
 
 
@@ -1104,4 +1105,5 @@
 
         self._order_ = self.default_order
+        print '.evolve: 2'
 
 
@@ -1127,20 +1129,25 @@
         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,5 +1209,5 @@
                 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,6 +1575,4 @@
         """
 
-        from Numeric import ones, sum, equal, Float
-
         N = len(self) # Number_of_triangles
         d = len(self.conserved_quantities)
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/ermapper_grids.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/ermapper_grids.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/ermapper_grids.py	(revision 5892)
@@ -2,7 +2,8 @@
 
 # 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,9 +12,9 @@
     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,5 +23,5 @@
     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,5 +58,5 @@
 
     # 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,5 +84,5 @@
     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,5 +164,5 @@
     return header                      
 
-def write_ermapper_data(grid, ofile, data_format = Numeric.Float32):
+def write_ermapper_data(grid, ofile, data_format = numpy.float32):
 
 
@@ -181,5 +182,5 @@
         
     
-    # 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,5 +194,5 @@
 
 
-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,6 +200,6 @@
     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
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 5892)
@@ -1,5 +1,3 @@
-
-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,6 +88,8 @@
         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,5 +138,5 @@
 
         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,11 +146,15 @@
                       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,12 +191,12 @@
             #   - 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,6 +278,6 @@
         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,4 +315,5 @@
             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,11 +321,11 @@
             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,5 +354,6 @@
 
         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,5 +382,5 @@
             indices = range(M)
 
-        return take(self.triangles, indices)
+        return numpy.take(self.triangles, indices, axis=0)
     
 
@@ -409,5 +415,6 @@
         
         #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,8 +446,8 @@
         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,5 +459,5 @@
                 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,5 +536,5 @@
 
         # 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,6 +542,6 @@
 
         # 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,8 +589,8 @@
         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,6 +605,3 @@
         """
 
-        return sum(self.areas)
-
-        
-        
+        return numpy.sum(self.areas)
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py	(revision 5892)
@@ -7,4 +7,6 @@
 from anuga.fit_interpolate.interpolate import Modeltime_too_late
 from anuga.fit_interpolate.interpolate import Modeltime_too_early
+
+import numpy
 
 
@@ -67,6 +69,5 @@
             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,8 +97,6 @@
             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,5 +135,4 @@
     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,5 +204,4 @@
 
         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,5 +219,5 @@
 
         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,7 +239,7 @@
             
             # 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
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py	(revision 5892)
@@ -2,4 +2,6 @@
 mesh file formats
 """
+
+import numpy
 
 
@@ -73,5 +75,4 @@
 
     from anuga.config import epsilon
-    from Numeric import zeros, Float, Int
 
     delta1 = float(len1)/m
@@ -93,5 +94,5 @@
     index = Index(n,m)
 
-    points = zeros( (Np,2), Float)
+    points = numpy.zeros( (Np,2), numpy.float)
 
     for i in range(m+1):
@@ -105,5 +106,5 @@
 
 
-    elements = zeros( (Nt,3), Int)
+    elements = numpy.zeros( (Nt,3), numpy.int)
     boundary = {}
     nt = -1
@@ -149,5 +150,4 @@
 
     from anuga.config import epsilon
-    from Numeric import zeros, Float, Int
 
     delta1 = float(len1)/m
@@ -208,5 +208,4 @@
     """
 
-    from Numeric import array
     import math
 
@@ -510,5 +509,5 @@
     """
 
-    from Numeric import array
+#    from numpy.oldnumeric import array
     import math
 
@@ -592,5 +591,4 @@
     """
 
-    from Numeric import array
     import math
 
@@ -686,9 +684,6 @@
     """
 
-    from Numeric import array
     import math
-
     from anuga.config import epsilon
-
 
     deltax = lenx/float(m)
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 5892)
@@ -11,6 +11,5 @@
 from anuga.caching import cache
 from math import pi, sqrt
-from Numeric import array, allclose
-        
+import numpy        
 
 class Mesh(General_mesh):
@@ -79,8 +78,4 @@
         """
 
-
-
-        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
-
         General_mesh.__init__(self, coordinates, triangles,
                               number_of_full_nodes=\
@@ -98,12 +93,12 @@
 
         #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,5 +119,5 @@
 
             #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,14 +128,14 @@
 
                 #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,7 +145,7 @@
                 #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,5 +383,4 @@
             self.element_tag is defined
         """
-        from Numeric import array, Int
 
         if tagged_elements is None:
@@ -395,5 +389,5 @@
             #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,5 +457,4 @@
         """
         
-        from Numeric import allclose, sqrt, array, minimum, maximum
         from anuga.utilities.numerical_tools import angle, ensure_numeric     
 
@@ -477,5 +470,5 @@
         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,10 +480,9 @@
             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,5 +585,5 @@
                 # 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,6 +621,4 @@
         from anuga.utilities.numerical_tools import anglediff
 
-        from Numeric import sort, allclose
-
         N = len(self)
 
@@ -672,9 +662,9 @@
 
                 # 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,6 +720,6 @@
 
         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,5 +732,5 @@
                 count[i] += 1
 
-        assert allclose(count, self.number_of_triangles_per_node)
+        assert numpy.allclose(count, self.number_of_triangles_per_node)
 
 
@@ -805,5 +795,4 @@
         """
 
-        from Numeric import arange
         from anuga.utilities.numerical_tools import histogram, create_bins
 
@@ -1141,8 +1130,8 @@
 
             # 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,7 +1146,7 @@
             # 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,5 +1228,5 @@
         assert isinstance(segment, Triangle_intersection), msg
         
-        midpoint = sum(array(segment.segment))/2
+        midpoint = numpy.sum(numpy.array(segment.segment))/2
         midpoints.append(midpoint)
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/pmesh2domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/pmesh2domain.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/pmesh2domain.py	(revision 5892)
@@ -9,4 +9,5 @@
 import sys
 
+import numpy
 
 
@@ -161,5 +162,4 @@
     """
 
-    from Numeric import transpose
     from load_mesh.loadASCII import import_mesh_file
 
@@ -172,9 +172,10 @@
     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)
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5892)
@@ -15,6 +15,5 @@
 """
 
-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,7 +37,7 @@
         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,13 +56,13 @@
 
         # 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,5 +71,5 @@
         # 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,7 +77,7 @@
 
         # 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,5 +381,4 @@
         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,5 +446,5 @@
 
         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,5 +456,5 @@
                 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,6 +472,6 @@
                                                      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,10 +608,8 @@
         """
 
-        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,8 +636,8 @@
 
         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,9 +672,8 @@
         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,5 +711,6 @@
                 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,6 +776,6 @@
 
 
-        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,7 +909,7 @@
         # 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,6 +1114,4 @@
         """
         
-        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,7 +1137,7 @@
             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,9 +1145,9 @@
             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,9 +1172,9 @@
                     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,8 +1190,6 @@
         """
 
-        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,7 +1265,4 @@
         """
 
-        from Numeric import concatenate, zeros, Float, Int, array, reshape
-
-
         if smooth is None:
             # Take default from domain
@@ -1284,5 +1273,5 @@
 
         if precision is None:
-            precision = Float
+            precision = numpy.float
             
 
@@ -1293,5 +1282,5 @@
             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,5 +1332,5 @@
             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)
 
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 5892)
@@ -11,5 +11,5 @@
 
 #include "Python.h"
-#include "Numeric/arrayobject.h"
+#include "numpy/arrayobject.h"
 #include "math.h"
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/region.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/region.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/region.py	(revision 5892)
@@ -8,5 +8,7 @@
 # FIXME (DSG-DSG) add better comments
 
-from Numeric import average
+import numpy
+
+
 class Region:
     """Base class for modifying quantities based on a region.
@@ -101,7 +103,7 @@
                 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:
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py	(revision 5892)
@@ -18,6 +18,5 @@
 
 from mesh_factory import rectangular
-from Numeric import array
-    
+
 
 ######################
@@ -78,6 +77,4 @@
 domain.set_quantity('stage', Z)
 
-from Numeric import allclose
-
 #Evolve
 for t in domain.evolve(yieldstep = 0.1, finaltime = 30):
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py	(revision 5892)
@@ -6,5 +6,5 @@
 from domain import *
 from anuga.config import epsilon
-from Numeric import allclose, array, ones, Float
+import numpy
 
 
@@ -64,5 +64,5 @@
 
 
-        assert domain.get_conserved_quantities(0, edge=1) == 0.
+        assert numpy.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
 
 
@@ -95,43 +95,43 @@
         #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,6 +179,6 @@
         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,5 +188,5 @@
         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,6 +314,6 @@
         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,5 +382,5 @@
         domain.check_integrity()
 
-        assert allclose(domain.neighbours, [[-1,-2,-3]])
+        assert numpy.allclose(domain.neighbours, [[-1,-2,-3]])
 
 
@@ -471,6 +471,6 @@
                                       [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,9 +484,9 @@
 
         #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,6 +527,6 @@
 
         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,16 +535,16 @@
         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,6 +578,6 @@
                                       [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,9 +591,9 @@
 
         #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,5 +610,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -627,5 +626,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -645,21 +643,21 @@
         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,5 +668,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -690,31 +687,31 @@
         
         #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]])
 
 #-------------------------------------------------------------
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_ermapper.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_ermapper.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_ermapper.py	(revision 5892)
@@ -5,5 +5,5 @@
 
 import ermapper_grids
-import Numeric
+import numpy
 from os import remove
 
@@ -19,5 +19,5 @@
         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,5 +25,5 @@
         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,5 +31,5 @@
         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,14 +40,14 @@
         # 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,5 +57,5 @@
         # 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,5 +66,5 @@
 
         # 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,5 +75,5 @@
 
         # 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)
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 5892)
@@ -7,5 +7,5 @@
 
 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,5 +23,4 @@
     def test_get_vertex_coordinates(self):
         from mesh_factory import rectangular
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -30,5 +29,5 @@
 
 
-        assert allclose(domain.get_nodes(), nodes)
+        assert numpy.allclose(domain.get_nodes(), nodes)
 
 
@@ -41,5 +40,5 @@
             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,27 +54,27 @@
         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,5 +85,4 @@
         """
         from mesh_factory import rectangular
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -93,5 +91,5 @@
 
 
-        assert allclose(domain.get_nodes(), nodes)
+        assert numpy.allclose(domain.get_nodes(), nodes)
 
 
@@ -104,5 +102,5 @@
             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,5 +112,4 @@
         """
         from mesh_factory import rectangular
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -121,7 +118,7 @@
 
         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,16 +127,15 @@
         """
         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,5 +158,5 @@
             #print count
             #
-            assert allclose(count, domain.number_of_triangles_per_node)
+            assert numpy.allclose(count, domain.number_of_triangles_per_node)
             
             # Check indices
@@ -196,7 +192,7 @@
         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,7 +200,8 @@
         # 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,5 +210,5 @@
         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,5 +221,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -239,5 +235,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -273,21 +268,21 @@
         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,19 +305,19 @@
         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,5 +328,4 @@
 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)
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py	(revision 5892)
@@ -6,5 +6,5 @@
 from generic_boundary_conditions import *
 from anuga.config import epsilon
-from Numeric import allclose, array
+import numpy
 
 
@@ -44,5 +44,5 @@
 
         q = Bd.evaluate()
-        assert allclose(q, x)
+        assert numpy.allclose(q, x)
 
 
@@ -98,5 +98,5 @@
         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,14 +175,14 @@
 
         #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,10 +193,10 @@
         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])
 
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_ghost.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_ghost.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_ghost.py	(revision 5892)
@@ -6,7 +6,5 @@
 from domain import *
 from anuga.config import epsilon
-from Numeric import allclose, array, ones, Float
-
-
+import numpy
 
 
@@ -43,5 +41,5 @@
 
 
-        assert domain.get_conserved_quantities(0, edge=1) == 0.
+        assert numpy.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
 
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py	(revision 5892)
@@ -1,5 +1,3 @@
 #!/usr/bin/env python
-
-
 
 #FIXME: Seperate the tests for mesh and general_mesh
@@ -13,5 +11,5 @@
 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,5 +18,5 @@
 
 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,14 +61,14 @@
         #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,16 +79,16 @@
 
         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,8 +235,8 @@
 
         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,8 +250,8 @@
 
         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,5 +267,5 @@
         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,8 +301,8 @@
         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,6 +695,6 @@
         get values based on triangle lists.
         """
+
         from mesh_factory import rectangular
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -716,6 +714,4 @@
     def test_boundary_polygon(self):
         from mesh_factory import rectangular
-        #from mesh import Mesh
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -727,7 +723,7 @@
 
         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,7 +732,4 @@
 
     def test_boundary_polygon_II(self):
-        from Numeric import zeros, Float
-        
-
         #Points
         a = [0.0, 0.0] #0
@@ -764,5 +757,5 @@
 
         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,7 +767,4 @@
         """Same as II but vertices ordered differently
         """
-
-        from Numeric import zeros, Float
-
 
         #Points
@@ -804,5 +794,5 @@
 
         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,7 +805,4 @@
 	is partitioned using pymetis.
         """
-
-        from Numeric import zeros, Float
-
 
         #Points
@@ -847,5 +834,5 @@
 	
 	# 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,5 +851,4 @@
         """
 
-        from Numeric import zeros, Float
         from mesh_factory import rectangular        
 
@@ -907,6 +893,4 @@
         
         """
-        from Numeric import zeros, Float
-        
 
         #Points
@@ -955,6 +939,6 @@
         
         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,5 +982,5 @@
                   [  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,5 +1085,5 @@
                   [  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,5 +1126,5 @@
             assert is_inside_polygon(p, P)            
 
-        assert allclose(P, Pref)    
+        assert numpy.allclose(P, Pref)    
 
     def test_lone_vertices(self):
@@ -1194,5 +1178,5 @@
         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,5 +1225,5 @@
 
         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,14 +1274,14 @@
             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,5 +1290,5 @@
 
             msg = 'Segments do not add up'
-            assert allclose(total_length, 2), msg
+            assert numpy.allclose(total_length, 2), msg
             
 
@@ -1341,10 +1325,10 @@
         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,5 +1339,5 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 2), msg
+        assert numpy.allclose(total_length, 2), msg
         
 
@@ -1395,10 +1379,10 @@
         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,5 +1392,5 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 1.5), msg            
+        assert numpy.allclose(total_length, 1.5), msg            
 
 
@@ -1442,7 +1426,7 @@
         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,5 +1435,5 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 4*s2), msg
+        assert numpy.allclose(total_length, 4*s2), msg
 
 
@@ -1466,7 +1450,7 @@
         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,5 +1459,5 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 4*s2), msg                    
+        assert numpy.allclose(total_length, 4*s2), msg                    
 
 
@@ -1491,7 +1475,7 @@
         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,5 +1484,5 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 4*s2), msg                        
+        assert numpy.allclose(total_length, 4*s2), msg                        
 
 
@@ -1515,7 +1499,7 @@
         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,5 +1508,5 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 4*s2), msg                        
+        assert numpy.allclose(total_length, 4*s2), msg                        
 
 
@@ -1540,7 +1524,7 @@
         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,5 +1533,5 @@
 
         msg = 'Segments do not add up'
-        assert allclose(total_length, 2*s2), msg
+        assert numpy.allclose(total_length, 2*s2), msg
 
 
@@ -1562,6 +1546,6 @@
         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,6 +1582,6 @@
         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,6 +1594,6 @@
         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,6 +1607,6 @@
         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,10 +1642,10 @@
         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,18 +1680,18 @@
         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,25 +1731,25 @@
         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,6 +1759,6 @@
                     
             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,5 +1826,5 @@
             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)
 
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py	(revision 5892)
@@ -3,7 +3,5 @@
 
 import unittest
-
-from Numeric import allclose, array
-
+import numpy
 
 #from anuga.pyvolution.pmesh2domain import *
@@ -66,7 +64,7 @@
 
          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,22 +78,22 @@
          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,7 +157,7 @@
         
          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,22 +172,22 @@
          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,7 +226,7 @@
          #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 = {}
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 5892)
@@ -7,5 +7,5 @@
 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,5 +18,5 @@
 #Aux for fit_interpolate.fit example
 def linear_function(point):
-    point = array(point)
+    point = numpy.array(point)
     return point[:,0]+point[:,1]
 
@@ -63,5 +63,5 @@
 
         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,9 +84,9 @@
 
         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,7 +94,7 @@
     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,5 +102,5 @@
         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,25 +108,25 @@
 
         #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,10 +148,10 @@
 
         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,14 +192,14 @@
 
         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,9 +212,9 @@
 
         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,5 +227,5 @@
 
         v = quantity.get_values(interpolation_points = [[x,y]])
-        assert allclose(v, 6)        
+        assert numpy.allclose(v, 6)        
 
         
@@ -244,26 +244,26 @@
         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,16 +288,16 @@
 
         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,15 +310,15 @@
         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,5 +331,5 @@
         #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,5 +342,5 @@
         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,12 +350,12 @@
         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,6 +365,6 @@
         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,14 +376,14 @@
 
 
-        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,6 +399,6 @@
 
         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,7 +408,7 @@
 
         # 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,7 +419,7 @@
         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,9 +429,9 @@
         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,9 +441,9 @@
         #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,12 +478,12 @@
 
         # 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,11 +491,11 @@
         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,9 +503,9 @@
         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,5 +540,5 @@
         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,5 +553,5 @@
         #print vertex_attributes
         quantity.set_values(vertex_attributes)
-        assert allclose(quantity.vertex_values.flat, answer)
+        assert numpy.allclose(quantity.vertex_values.ravel(), answer)
 
 
@@ -593,5 +593,5 @@
                           alpha = 0)
 
-        assert allclose( ref, [0,5,5] )
+        assert numpy.allclose( ref, [0,5,5] )
 
 
@@ -610,5 +610,5 @@
         #                    data_georef = data_georef,
         #                    alpha = 0)
-        assert allclose(quantity.vertex_values.flat, ref)
+        assert numpy.allclose(quantity.vertex_values.ravel(), ref)
 
 
@@ -621,5 +621,5 @@
 
         quantity.set_values(geospatial_data = geo, alpha = 0)
-        assert allclose(quantity.vertex_values.flat, ref)
+        assert numpy.allclose(quantity.vertex_values.ravel(), ref)
 
 
@@ -666,14 +666,14 @@
         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,5 +730,5 @@
         #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,16 +748,16 @@
         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,5 +815,5 @@
                             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,5 +821,5 @@
         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,13 +866,13 @@
         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,13 +910,13 @@
         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,7 +957,7 @@
                                       '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,12 +965,12 @@
                             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,7 +1039,7 @@
                                       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,12 +1047,12 @@
                             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,10 +1128,10 @@
         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,10 +1207,10 @@
 
 
-        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,26 +1226,26 @@
         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,11 +1262,11 @@
         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,7 +1276,7 @@
                              
         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,6 +1287,6 @@
         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,43 +1304,43 @@
         # 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,59 +1351,59 @@
         #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,6 +1431,6 @@
         #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,7 +1438,7 @@
         #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,7 +1446,7 @@
         #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,7 +1454,7 @@
         #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,11 +1464,11 @@
 
         #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,6 +1489,6 @@
         #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,7 +1496,7 @@
         #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,7 +1504,7 @@
         #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,7 +1512,7 @@
         #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,6 +1532,6 @@
         #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,7 +1540,7 @@
 
         #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,6 +1550,6 @@
 
         #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,5 +1558,5 @@
 
         #print quantity.vertex_values
-        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
+        assert numpy.allclose(quantity.centroid_values, quantity.centroid_backup_values)
 
 
@@ -1574,5 +1574,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
 
         #Extrapolate
@@ -1581,10 +1581,10 @@
         #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,8 +1614,7 @@
 
         #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,8 +1645,7 @@
 
         #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,8 +1674,7 @@
 
         #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,8 +1702,7 @@
 
         #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,5 +1714,5 @@
         #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,5 +1720,5 @@
         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,14 +1727,13 @@
         # 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,5 +1746,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
 
 
@@ -1760,8 +1755,8 @@
         #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,12 +1764,12 @@
         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,13 +1777,13 @@
         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,5 +1794,5 @@
         #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,5 +1800,5 @@
         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,14 +1808,14 @@
         #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,8 +1824,8 @@
         #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,9 +1833,9 @@
         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,11 +1845,11 @@
         #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,12 +1857,12 @@
         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,5 +1874,4 @@
         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,5 +1900,5 @@
 
         bed = domain.quantities['elevation'].vertex_values
-        stage = zeros(bed.shape, Float)
+        stage = numpy.zeros(bed.shape, numpy.float)
 
         h = 0.03
@@ -1929,36 +1923,36 @@
 
         #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,5 +1963,4 @@
         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,5 +1984,5 @@
 
         bed = domain.quantities['elevation'].vertex_values
-        stage = zeros(bed.shape, Float)
+        stage = numpy.zeros(bed.shape, numpy.float)
 
         h = 0.03
@@ -2005,8 +1998,8 @@
         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,14 +2014,14 @@
 
 
-        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,5 +2031,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2054,11 +2046,11 @@
         #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,5 +2060,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2087,5 +2078,5 @@
                             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,10 +2086,10 @@
                             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,9 +2098,9 @@
         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,7 +2112,7 @@
         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,5 +2121,5 @@
         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,10 +2133,10 @@
         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,5 +2146,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2171,7 +2161,7 @@
                             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,5 +2172,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2205,12 +2194,12 @@
 
         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,24 +2230,24 @@
         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,5 +2269,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2298,5 +2286,5 @@
         
         #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,6 +2299,6 @@
         #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,7 +2306,7 @@
         #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,5 +2333,5 @@
         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,5 +2339,5 @@
         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,5 +2368,5 @@
         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,5 +2375,5 @@
         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,10 +2381,10 @@
                [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,5 +2393,5 @@
                               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,5 +2404,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2438,6 +2425,6 @@
         #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,13 +2432,13 @@
         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,5 +2447,5 @@
         #print "subset",subset
         #print "answer",answer
-        assert allclose(subset, answer)
+        assert numpy.allclose(subset, answer)
 
         subset = quantity.get_values( indices=[1,5])
@@ -2466,5 +2453,5 @@
         #print "subset",subset
         #print "answer",answer
-        assert allclose(subset, answer)
+        assert numpy.allclose(subset, answer)
 
     def test_smooth_vertex_values(self):
@@ -2474,5 +2461,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -2501,6 +2487,6 @@
         
         #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,6 +2498,6 @@
                                 [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,10 +2508,6 @@
 #-------------------------------------------------------------
 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)
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py	(revision 5892)
@@ -6,6 +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,5 +44,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -64,33 +63,33 @@
         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,11 +97,11 @@
         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,11 +114,11 @@
         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,5 +128,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -147,11 +145,11 @@
         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,5 +160,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -180,11 +177,11 @@
 
         #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,5 +191,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -218,8 +214,8 @@
         #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,5 +225,4 @@
         from mesh_factory import rectangular
         from shallow_water import Domain
-        from Numeric import zeros, Float
 
         #Create basic mesh
@@ -252,12 +247,12 @@
         #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,4 +260,5 @@
 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)
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py	(revision 5892)
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 
-
 import unittest
-from Numeric import zeros, array, allclose, Float
+import numpy
+
 from math import sqrt, pi
 import tempfile, os
@@ -91,8 +91,8 @@
 
             #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,12 +100,12 @@
         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,5 +125,4 @@
         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,13 +147,19 @@
 
         # 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,4 +169,5 @@
              
             #domain1.write_time()
+        print 'test_util: 5'
 
 
@@ -182,8 +188,8 @@
 
         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,8 +200,8 @@
 
         #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,7 +215,7 @@
         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,8 +224,8 @@
 
         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,7 +237,7 @@
         #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,8 +246,8 @@
 
         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,7 +257,7 @@
         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,8 +267,8 @@
 
         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,8 +280,8 @@
         #
         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,7 +296,7 @@
         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,7 +310,7 @@
 
         #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,5 +332,4 @@
         from shallow_water import Domain, Dirichlet_boundary
         from mesh_factory import rectangular
-        from Numeric import take, concatenate, reshape
 
 
@@ -386,8 +391,8 @@
 
         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,8 +403,8 @@
 
         #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,7 +420,7 @@
 
         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,8 +429,8 @@
 
         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,7 +442,7 @@
         #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,8 +451,8 @@
 
         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,7 +462,7 @@
         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,8 +472,8 @@
 
         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,8 +485,8 @@
         #
         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,7 +501,7 @@
         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,7 +515,7 @@
 
         #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,5 +547,4 @@
         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,5 +588,5 @@
             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,5 +608,5 @@
 
         #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,5 +615,5 @@
                           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,5 +649,5 @@
                      self.failUnless( q == actual, 'Fail!')
                 else:
-                    assert allclose(q, actual)
+                    assert numpy.allclose(q, actual)
 
 
@@ -655,9 +659,9 @@
             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,5 +674,5 @@
                           quantities = domain.conserved_quantities,
                           interpolation_points = interpolation_points)
-        assert allclose(domain.starttime, start+delta)
+        assert numpy.allclose(domain.starttime, start+delta)
 
 
@@ -689,5 +693,5 @@
 
                 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,5 +707,4 @@
         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,5 +754,5 @@
             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,5 +777,5 @@
 
         #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,5 +784,5 @@
                           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,5 +820,5 @@
                      self.failUnless( q == actual, 'Fail!')
                 else:
-                    assert allclose(q, actual)
+                    assert numpy.allclose(q, actual)
 
         # now lets check points inside the mesh
@@ -863,5 +866,5 @@
                      self.failUnless( q == actual, 'Fail!')
                 else:
-                    assert allclose(q, actual)
+                    assert numpy.allclose(q, actual)
 
 
@@ -873,9 +876,9 @@
             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,5 +891,5 @@
                           quantities = domain.conserved_quantities,
                           interpolation_points = interpolation_points)
-        assert allclose(domain.starttime, start+delta)
+        assert numpy.allclose(domain.starttime, start+delta)
 
 
@@ -907,5 +910,5 @@
 
                 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,5 +958,5 @@
                           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,5 +965,5 @@
                           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,5 +972,5 @@
                           domain,
                           quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
-        assert allclose(domain.starttime, start+1)
+        assert numpy.allclose(domain.starttime, start+1)
 
         domain.starttime = start
@@ -987,8 +990,8 @@
 
             #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,12 +999,12 @@
         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,5 +1055,5 @@
         F = file_function(filename + '.tms', domain,
                           quantities = ['Attribute0', 'Attribute1', 'Attribute2'])        
-        assert allclose(domain.starttime, start+delta)
+        assert numpy.allclose(domain.starttime, start+delta)
 
 
@@ -1063,8 +1066,8 @@
 
             #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,12 +1075,12 @@
         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,24 +1094,22 @@
         #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,6 +1269,6 @@
         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,6 +1285,6 @@
         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,7 +1292,7 @@
         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,6 +1300,6 @@
         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,6 +1308,6 @@
         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,5 +1497,5 @@
             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,5 +1510,5 @@
             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,5 +1616,5 @@
             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,5 +1629,5 @@
             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,5 +1735,5 @@
             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,5 +1748,5 @@
             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,5 +1831,5 @@
 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)
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 5891)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 5892)
@@ -15,5 +15,5 @@
 
 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,5 +237,5 @@
     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,5 +303,5 @@
     # Get variables
     # if verbose: print 'Get variables'    
-    time = fid.variables['time'][:]    
+    time = numpy.array(fid.variables['time'][:]    )
 
     # Get time independent stuff
@@ -317,7 +317,7 @@
             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,5 +331,5 @@
             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,5 +351,5 @@
                 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,4 +373,8 @@
         # 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,5 +402,5 @@
         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,5 +524,5 @@
                 raise 'Illegal input to get_textual_float:', value
         else:
-            return format %float(value)
+            return format % float(value)
 
 
@@ -1036,5 +1040,5 @@
     """
 #    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,16 +1081,16 @@
     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,7 +1104,7 @@
     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,5 +1206,5 @@
                 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,5 +1505,5 @@
     
     # 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,5 +1539,5 @@
         #print "loners", loners
         #print "triangles before", triangles
-        triangles = choose(triangles,loners)
+        triangles = numpy.choose(triangles,loners)
         #print "triangles after", triangles
     return verts, triangles
@@ -1550,5 +1554,5 @@
 
         
-    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,5 +1836,5 @@
                 #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,5 +1947,5 @@
             #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,5 +2128,5 @@
     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,5 +2169,5 @@
         
     #convert to array for file_function
-    points_array = array(points,Float)
+    points_array = numpy.array(points, numpy.float)
         
     points_array = ensure_absolute(points_array)
