Changeset 5897

Timestamp:

Nov 6, 2008, 12:17:15 PM (15 years ago)

Author:

ole

Message:

Reverted numpy changes to the trunk that should have been made to the branch.
The command was svn merge -r 5895:5890 .

Location:

anuga_core

Files:

• compile_all.py (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/domain.py (modified) (15 diffs)
• source/anuga/abstract_2d_finite_volumes/ermapper_grids.py (modified) (9 diffs)
• source/anuga/abstract_2d_finite_volumes/general_mesh.py (modified) (16 diffs)
• source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (7 diffs)
• source/anuga/abstract_2d_finite_volumes/mesh_factory.py (modified) (9 diffs)
• source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py (modified) (20 diffs)
• source/anuga/abstract_2d_finite_volumes/pmesh2domain.py (modified) (3 diffs)
• source/anuga/abstract_2d_finite_volumes/quantity.py (modified) (24 diffs)
• source/anuga/abstract_2d_finite_volumes/quantity_ext.c (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/region.py (modified) (2 diffs)
• source/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py (modified) (2 diffs)
• source/anuga/abstract_2d_finite_volumes/test_domain.py (modified) (18 diffs)
• source/anuga/abstract_2d_finite_volumes/test_ermapper.py (modified) (8 diffs)
• source/anuga/abstract_2d_finite_volumes/test_general_mesh.py (modified) (20 diffs)
• source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py (modified) (5 diffs)
• source/anuga/abstract_2d_finite_volumes/test_ghost.py (modified) (2 diffs)
• source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py (modified) (51 diffs)
• source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py (modified) (6 diffs)
• source/anuga/abstract_2d_finite_volumes/test_quantity.py (modified) (124 diffs)
• source/anuga/abstract_2d_finite_volumes/test_region.py (modified) (14 diffs)
• source/anuga/abstract_2d_finite_volumes/test_util.py (modified) (67 diffs)
• source/anuga/abstract_2d_finite_volumes/util.py (modified) (20 diffs)
• source/anuga/caching/caching.py (modified) (6 diffs)
• source/anuga/caching/test_caching.py (modified) (7 diffs)
• source/anuga/utilities/cg_solve.py (modified) (5 diffs)
• source/anuga/utilities/compile.py (modified) (2 diffs)
• source/anuga/utilities/interp.py (modified) (5 diffs)
• source/anuga/utilities/numerical_tools.py (modified) (9 diffs)
• source/anuga/utilities/polygon.py (modified) (25 diffs)
• source/anuga/utilities/polygon_ext.c (modified) (4 diffs)
• source/anuga/utilities/quad.py (modified) (4 diffs)
• source/anuga/utilities/sparse.py (modified) (15 diffs)
• source/anuga/utilities/sparse_ext.c (modified) (1 diff)
• source/anuga/utilities/test_cg_solve.py (modified) (12 diffs)
• source/anuga/utilities/test_data_audit.py (modified) (1 diff)
• source/anuga/utilities/test_numerical_tools.py (modified) (8 diffs)
• source/anuga/utilities/test_polygon.py (modified) (41 diffs)
• source/anuga/utilities/test_quad.py (modified) (3 diffs)
• source/anuga/utilities/test_sparse.py (modified) (11 diffs)
• source/anuga/utilities/test_system_tools.py (modified) (4 diffs)
• source/anuga/utilities/test_xml_tools.py (modified) (1 diff)
• source/anuga/utilities/util_ext.c (modified) (1 diff)
• source/anuga/utilities/util_ext.h (modified) (1 diff)

anuga_core/compile_all.py

@@ -3 +3 @@
 buildroot = os.getcwd()
 
-os.chdir('source')
+os.chdir('source_numpy_conversion')
 
 os.chdir('anuga')

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -8 +8 @@
 """
 
-##from numpy.oldnumeric import allclose, argmax, zeros, Float
-import numpy
-
+from Numeric import allclose, argmax, zeros, Float
 from anuga.config import epsilon
 from anuga.config import beta_euler, beta_rk2
@@ -107 +105 @@
 
         if verbose: print 'Initialising Domain'
+        from Numeric import zeros, Float, Int, ones
         from quantity import Quantity
 
@@ -156 +155 @@
         for key in self.full_send_dict:
             buffer_shape = self.full_send_dict[key][0].shape[0]
-            self.full_send_dict[key].append(numpy.zeros( (buffer_shape,self.nsys), numpy.float))
+            self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
 
 
         for key in self.ghost_recv_dict:
             buffer_shape = self.ghost_recv_dict[key][0].shape[0]
-            self.ghost_recv_dict[key].append(numpy.zeros( (buffer_shape,self.nsys), numpy.float))
+            self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
 
 
@@ -169 +168 @@
         N = len(self) #number_of_elements
         self.number_of_elements = N
-        self.tri_full_flag = numpy.ones(N, numpy.int)
+        self.tri_full_flag = ones(N, Int)
         for i in self.ghost_recv_dict.keys():
             for id in self.ghost_recv_dict[i][0]:
@@ -176 +175 @@
         # Test the assumption that all full triangles are store before
         # the ghost triangles.
-        if not numpy.allclose(self.tri_full_flag[:self.number_of_full_nodes],1):
+        if not allclose(self.tri_full_flag[:self.number_of_full_nodes],1):
             print 'WARNING:  Not all full triangles are store before ghost triangles'
                         
@@ -231 +230 @@
         # calculation
         N = len(self) # Number_of_triangles
-        self.already_computed_flux = numpy.zeros((N, 3), numpy.int)
+        self.already_computed_flux = zeros((N, 3), 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 = numpy.zeros(N, numpy.float)
+        self.max_speed = zeros(N, Float)
 
         if mesh_filename is not None:
@@ -267 +266 @@
         """
 
+        from Numeric import zeros, Float
+
         if not (vertex is None or edge is None):
             msg = 'Values for both vertex and edge was specified.'
@@ -272 +273 @@
             raise msg
 
-        q = numpy.zeros( len(self.conserved_quantities), numpy.float)
+        q = zeros( len(self.conserved_quantities), Float)
 
         for i, name in enumerate(self.conserved_quantities):
@@ -768 +769 @@
             # Find index of largest computed flux speed
             if triangle_id is None:
-                k = self.k = numpy.argmax(self.max_speed)
+                k = self.k = argmax(self.max_speed)
             else:
                 errmsg = 'Triangle_id %d does not exist in mesh: %s' %(triangle_id,
@@ -1085 +1086 @@
 
         """
-        print '.evolve: 0'
 
         from anuga.config import min_timestep, max_timestep, epsilon
@@ -1096 +1096 @@
                %self.get_boundary_tags()
         assert hasattr(self, 'boundary_objects'), msg
-        print '.evolve: 1'
 
 
@@ -1105 +1104 @@
 
         self._order_ = self.default_order
-        print '.evolve: 2'
 
 
@@ -1129 +1127 @@
         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:
@@ -1209 +1202 @@
                 self.number_of_steps = 0
                 self.number_of_first_order_steps = 0
-                self.max_speed = numpy.zeros(N, numpy.float)
+                self.max_speed = zeros(N, Float)
 
     def evolve_one_euler_step(self, yieldstep, finaltime):
@@ -1575 +1568 @@
         """
 
+        from Numeric import ones, sum, equal, Float
+
         N = len(self) # Number_of_triangles
         d = len(self.conserved_quantities)

anuga_core/source/anuga/abstract_2d_finite_volumes/ermapper_grids.py

@@ -2 +2 @@
 
 # from os import open, write, read
-##import numpy.oldnumeric as Numeric
-import numpy
-
-celltype_map = {'IEEE4ByteReal': numpy.float32, 'IEEE8ByteReal': numpy.float64}
+import Numeric
+
+celltype_map = {'IEEE4ByteReal': Numeric.Float32, 'IEEE8ByteReal': Numeric.Float64}
 
 
@@ -12 +11 @@
     write_ermapper_grid(ofile, data, header = {}):
     
-    Function to write a 2D NumPy array to an ERMapper grid.  There are a series of conventions adopted within
+    Function to write a 2D Numeric 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 NumPy 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 Numeric 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
@@ -23 +22 @@
     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:       NumPy.array - 2D array containing the data to be output to the grid
+    data:       Numeric.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"')
@@ -58 +57 @@
 
     # Check that the data is a 2 dimensional array
-    data_size = numpy.shape(data)
+    data_size = Numeric.shape(data)
     assert len(data_size) == 2
     
@@ -84 +83 @@
     nrofcellsperlines = int(header['nrofcellsperline'])
     data = read_ermapper_data(data_file)
-    data = numpy.reshape(data,(nroflines,nrofcellsperlines))
+    data = Numeric.reshape(data,(nroflines,nrofcellsperlines))
     return data
     
@@ -164 +163 @@
     return header                      
 
-def write_ermapper_data(grid, ofile, data_format = numpy.float32):
+def write_ermapper_data(grid, ofile, data_format = Numeric.Float32):
 
 
@@ -182 +181 @@
         
     
-    # 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)
 
@@ -194 +193 @@
 
 
-def read_ermapper_data(ifile, data_format = numpy.float32):
+def read_ermapper_data(ifile, data_format = Numeric.Float32):
     # open input file in a binary format and read the input string
     fid = open(ifile,'rb')
@@ -200 +199 @@
     fid.close()
 
-    # convert input string to required format (Note default format is numpy.float32)
-    grid_as_float = numpy.fromstring(input_string,data_format)
+    # convert input string to required format (Note default format is Numeric.Float32)
+    grid_as_float = Numeric.fromstring(input_string,data_format)
     return grid_as_float
 

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -1 @@
-import numpy
+
+from Numeric import concatenate, reshape, take, allclose
+from Numeric import array, zeros, Int, Float, sqrt, sum, arange
 
 from anuga.coordinate_transforms.geo_reference import Geo_reference
@@ -88 +90 @@
         if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' 
 
-##        self.triangles = array(triangles, Int)
-##        self.nodes = array(nodes, Float)
-        self.triangles = numpy.array(triangles, numpy.int)
-        self.nodes = numpy.array(nodes, numpy.float)
+        self.triangles = array(triangles, Int)
+        self.nodes = array(nodes, Float)
 
 
@@ -138 +138 @@
 
         msg = 'Vertex indices reference non-existing coordinate sets'
-        assert max(self.triangles.ravel()) < self.nodes.shape[0], msg
+        assert max(self.triangles.flat) < self.nodes.shape[0], msg
 
 
@@ -146 +146 @@
                       max(self.nodes[:,0]), max(self.nodes[:,1]) ]
 
-##        self.xy_extent = array(xy_extent, Float)
-        self.xy_extent = numpy.array(xy_extent, numpy.float)
+        self.xy_extent = array(xy_extent, 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 = numpy.zeros((N, 6), numpy.float)
-        self.areas = numpy.zeros(N, numpy.float)
-        self.edgelengths = numpy.zeros((N, 3), numpy.float)
+        self.normals = zeros((N, 6), Float)
+        self.areas = zeros(N, Float)
+        self.edgelengths = zeros((N, 3), Float)
 
         # Get x,y coordinates for all triangles and store
@@ -191 +187 @@
             #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
 
-            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))
+            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))
 
             # Normalise
@@ -278 +274 @@
         if absolute is True:
             if not self.geo_reference.is_absolute():
-                return V + numpy.array([self.geo_reference.get_xllcorner(),
-                                        self.geo_reference.get_yllcorner()])
+                return V + array([self.geo_reference.get_xllcorner(),
+                                  self.geo_reference.get_yllcorner()])
             else:
                 return V
@@ -320 +316 @@
             i3 = 3*i  
             if absolute is True and not self.geo_reference.is_absolute():
-                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])
+                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])
             else:
-                return numpy.array([V[i3,:], V[i3+1,:], V[i3+2,:]])
+                return array([V[i3,:], V[i3+1,:], V[i3+2,:]])
                 
 
@@ -353 +349 @@
 
         M = self.number_of_triangles
-##        vertex_coordinates = zeros((3*M, 2), Float)
-        vertex_coordinates = numpy.zeros((3*M, 2), numpy.float)
+        vertex_coordinates = zeros((3*M, 2), Float)
 
         for i in range(M):
@@ -381 +376 @@
             indices = range(M)
 
-        return numpy.take(self.triangles, indices, axis=0)
+        return take(self.triangles, indices)
     
 
@@ -414 +409 @@
         
         #T = reshape(array(range(K)).astype(Int), (M,3))
-##        T = reshape(arange(K).astype(Int), (M,3))  # Faster
-        T = numpy.reshape(numpy.arange(K).astype(numpy.int), (M,3))  # Faster
+        T = reshape(arange(K).astype(Int), (M,3))  # Faster
         
         return T     
@@ -445 +439 @@
         if node is not None:
             # Get index for this node
-            first = numpy.sum(self.number_of_triangles_per_node[:node])
+            first = sum(self.number_of_triangles_per_node[:node])
             
             # Get number of triangles for this node
-            count = int(self.number_of_triangles_per_node[node])
+            count = self.number_of_triangles_per_node[node]
 
             for i in range(count):
@@ -458 +452 @@
                 triangle_list.append( (volume_id, vertex_id) )
 
-            triangle_list = numpy.array(triangle_list)    
+            triangle_list = array(triangle_list)    
         else:
             # Get info for all nodes recursively.
@@ -535 +529 @@
 
         # Count number of triangles per node
-        number_of_triangles_per_node = numpy.zeros(self.number_of_full_nodes)
+        number_of_triangles_per_node = zeros(self.number_of_full_nodes)
         for volume_id, triangle in enumerate(self.get_triangles()):
             for vertex_id in triangle:
@@ -541 +535 @@
 
         # Allocate space for inverted structure
-        number_of_entries = numpy.sum(number_of_triangles_per_node)
-        vertex_value_indices = numpy.zeros(number_of_entries, numpy.int)
+        number_of_entries = sum(number_of_triangles_per_node)
+        vertex_value_indices = zeros(number_of_entries)
 
         # Register (triangle, vertex) indices for each node
@@ -588 +582 @@
         Y = C[:,1:6:2].copy()
 
-        xmin = min(X.ravel())
-        xmax = max(X.ravel())
-        ymin = min(Y.ravel())
-        ymax = max(Y.ravel())
+        xmin = min(X.flat)
+        xmax = max(X.flat)
+        ymin = min(Y.flat)
+        ymax = max(Y.flat)
         #print "C",C
         return xmin, xmax, ymin, ymax
@@ -604 +598 @@
         """
 
-        return numpy.sum(self.areas)
+        return sum(self.areas)
+
+        
+        

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -7 +7 @@
 from anuga.fit_interpolate.interpolate import Modeltime_too_late
 from anuga.fit_interpolate.interpolate import Modeltime_too_early
-
-import numpy
 
 
@@ -69 +67 @@
             raise Exception, msg
 
-        self.conserved_quantities=numpy.array(conserved_quantities).astype(numpy.float)
+        from Numeric import array, Float
+        self.conserved_quantities=array(conserved_quantities).astype(Float)
 
     def __repr__(self):
@@ -97 +96 @@
             raise msg
 
+
+        from Numeric import array, Float
         try:
-            q = numpy.array(q).astype(numpy.float)
+            q = array(q).astype(Float)
         except:
             msg = 'Return value from time boundary function could '
@@ -135 +136 @@
     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
@@ -204 +206 @@
 
         import time
+        from Numeric import array, zeros, Float
         from anuga.config import time_format
         from anuga.abstract_2d_finite_volumes.util import file_function
@@ -219 +222 @@
 
         if verbose: print 'Find midpoint coordinates of entire boundary'
-        self.midpoint_coordinates = numpy.zeros( (len(domain.boundary), 2), numpy.float)
+        self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
         boundary_keys = domain.boundary.keys()
 
@@ -239 +242 @@
             
             # Compute midpoints
-            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])
+            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])
 
             # Convert to absolute UTM coordinates

anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/pmesh2domain.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

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

anuga_core/source/anuga/abstract_2d_finite_volumes/region.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_ermapper.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_ghost.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py

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

anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

@@ -1 +1 @@
 #!/usr/bin/env python
 
+
 import unittest
-import numpy
-
+from Numeric import zeros, array, allclose, Float
 from math import sqrt, pi
 import tempfile, os
@@ -91 +91 @@
 
             #Exact linear intpolation
-            assert numpy.allclose(q[0], 2*t)
+            assert allclose(q[0], 2*t)
             if i%6 == 0:
-                assert numpy.allclose(q[1], t**2)
-                assert numpy.allclose(q[2], sin(t*pi/600))
+                assert allclose(q[1], t**2)
+                assert allclose(q[2], sin(t*pi/600))
 
         #Check non-exact
@@ -100 +100 @@
         t = 90 #Halfway between 60 and 120
         q = F(t)
-        assert numpy.allclose( (120**2 + 60**2)/2, q[1] )
-        assert numpy.allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
+        assert allclose( (120**2 + 60**2)/2, q[1] )
+        assert 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 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] )
+        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] )
 
         os.remove(filename + '.txt')
@@ -125 +125 @@
         from shallow_water import Domain, Dirichlet_boundary
         from mesh_factory import rectangular
+        from Numeric import take, concatenate, reshape
 
         #Create basic mesh and shallow water domain
@@ -147 +148 @@
 
         # 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:
@@ -169 +164 @@
              
             #domain1.write_time()
-        print 'test_util: 5'
 
 
@@ -188 +182 @@
 
         last_time_index = len(time)-1 #Last last_time_index
-        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)
+        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)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -200 +194 @@
 
         #And the midpoints are found now
-        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)
+        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)
         d_midpoints = (diag[1:] + diag[:-1])/2
 
@@ -215 +209 @@
         assert not T[-1] == T[-2], msg
         t = time[last_time_index]
-        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)
+        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)
 
 
@@ -224 +218 @@
 
         timestep = 0 #First timestep
-        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)
+        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)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -237 +231 @@
         #Let us see if the file function can find the correct
         #values
-        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)
+        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)
 
 
@@ -246 +240 @@
 
         timestep = 33
-        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)
+        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)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -257 +251 @@
         r2 = (D[2] + D[3])/2
 
-        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)
+        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)
 
 
@@ -267 +261 @@
 
         timestep = 15
-        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)
+        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)
 
         #Reference interpolated values at midpoints on diagonal at
@@ -280 +274 @@
         #
         timestep = 16
-        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