Changeset 3945

Timestamp:

Nov 8, 2006, 5:35:23 PM (18 years ago)

Author:

ole

Message:

One large step towards major cleanup. This has mainly to do with
the way vertex coordinates are handled internally.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/general_mesh.py (modified) (10 diffs)
• abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (1 diff)
• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (2 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (7 diffs)
• abstract_2d_finite_volumes/test_general_mesh.py (modified) (2 diffs)
• abstract_2d_finite_volumes/test_neighbour_mesh.py (modified) (3 diffs)
• abstract_2d_finite_volumes/test_quantity.py (modified) (5 diffs)
• fit_interpolate/fit.py (modified) (4 diffs)
• fit_interpolate/interpolate.py (modified) (1 diff)
• utilities/numerical_tools.py (modified) (1 diff)
• utilities/quad.py (modified) (3 diffs)
• utilities/test_quad.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -6 +6 @@
 
 class General_mesh:
-    """Collection of triangular elements (purely geometric)
+    """Collection of 2D triangular elements
 
     A triangular element is defined in terms of three vertex ids,
-    ordered counter clock-wise,
-    each corresponding to a given coordinate set.
-    Vertices from different elements can point to the same
-    coordinate set.
-
-    Coordinate sets are implemented as an N x 2 Numeric array containing
+    ordered counter clock-wise, each corresponding to a given node
+    which is represented as a coordinate set (x,y).
+    Vertices from different triangles can point to the same node.
+    The nodes are implemented as an Nx2 Numeric array containing the
     x and y coordinates.
 
 
     To instantiate:
-       Mesh(coordinates, triangles)
+       Mesh(nodes, triangles)
 
     where
 
-      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
+      nodes is either a list of 2-tuples or an Nx2 Numeric array of
       floats representing all x, y coordinates in the mesh.
 
-      triangles is either a list of 3-tuples or an Nx3 Numeric array of
+      triangles is either a list of 3-tuples or an Mx3 Numeric array of
       integers representing indices of all vertices in the mesh.
-      Each vertex is identified by its index i in [0, M-1].
+      Each vertex is identified by its index i in [0, N-1].
 
 
     Example:
+
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -37 +36 @@
         e = [2.0, 2.0]
 
-        points = [a, b, c, e]
-        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
-        mesh = Mesh(points, triangles)
-
-        #creates two triangles: bac and bce
+        nodes = [a, b, c, e]
+        triangles = [ [1,0,2], [1,2,3] ]   # bac, bce
+
+        # Create mesh with two triangles: bac and bce    
+        mesh = Mesh(nodes, triangles)
+
+
 
     Other:
 
-      In addition mesh computes an Nx6 array called vertex_coordinates.
+      In addition mesh computes an Mx6 array called vertex_coordinates.
       This structure is derived from coordinates and contains for each
       triangle the three x,y coordinates at the vertices.
 
-
-        This is a cut down version of mesh from mesh.py
+      See neighbourmesh.py for a specialisation of the general mesh class
+      which includes information about neighbours and the mesh boundary.
+
+      The mesh object is purely geometrical and contains no information
+      about quantities defined on the mesh.
+
     """
 
     #FIXME: It would be a good idea to use geospatial data as an alternative
     #input
-    def __init__(self, coordinates, triangles,
+    def __init__(self, nodes, triangles,
+                 geo_reference=None,                 
                  number_of_full_nodes=None,
                  number_of_full_triangles=None,                 
-                 geo_reference=None,
                  verbose=False):
-        """
-        Build triangles from x,y coordinates (sequence of 2-tuples or
-        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
-        or Nx3 Numeric array of non-negative integers).
-
-        origin is a 3-tuple consisting of UTM zone, easting and northing.
-        If specified coordinates are assumed to be relative to this origin.
+        """Build triangular 2d mesh from nodes and triangle information
+
+        Input:
+        
+          nodes: x,y coordinates represented as a sequence of 2-tuples or
+                 a Nx2 Numeric array of floats.
+                 
+          triangles: sequence of 3-tuples or Mx3 Numeric array of
+                     non-negative integers representing indices into
+                     the nodes array.
+       
+          georeference (optional): If specified coordinates are
+          assumed to be relative to this origin.
 
 
@@ -74 +85 @@
         In this case it is usefull to specify the number of real (called full)
         nodes and triangles. If omitted they will default to all.
+          
         """
 
         if verbose: print 'General_mesh: Building basic mesh structure' 
 
-        self.triangles = array(triangles,Int)
-        self.coordinates = array(coordinates,Float)
+        self.triangles = array(triangles, Int)
+        self.nodes = array(nodes, Float)
+
 
         # Register number of elements and nodes 
         self.number_of_triangles = N = self.triangles.shape[0]
-        self.number_of_nodes = self.coordinates.shape[0]        
+        self.number_of_nodes = self.nodes.shape[0]        
 
         
 
         if number_of_full_nodes is None:
-            self.number_of_full_nodes=self.number_of_nodes
+            self.number_of_full_nodes = self.number_of_nodes
         else:
             assert int(number_of_full_nodes)
-            self.number_of_full_nodes=number_of_full_nodes            
+            self.number_of_full_nodes = number_of_full_nodes            
 
 
         if number_of_full_triangles is None:
-            self.number_of_full_triangles=self.number_of_triangles           
+            self.number_of_full_triangles = self.number_of_triangles           
         else:
             assert int(number_of_full_triangles)            
-            self.number_of_full_triangles=number_of_full_triangles
+            self.number_of_full_triangles = number_of_full_triangles
         
         
@@ -114 +127 @@
 
         # Input checks
-        msg = 'Triangles must an Nx3 Numeric array or a sequence of 3-tuples. '
+        msg = 'Triangles must an Mx3 Numeric array or a sequence of 3-tuples. '
         msg += 'The supplied array has the shape: %s'\
                %str(self.triangles.shape)
         assert len(self.triangles.shape) == 2, msg
 
-        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
+        msg = 'Nodes must an Nx2 Numeric array or a sequence of 2-tuples'
         msg += 'The supplied array has the shape: %s'\
-               %str(self.coordinates.shape)
-        assert len(self.coordinates.shape) == 2, msg
+               %str(self.nodes.shape)
+        assert len(self.nodes.shape) == 2, msg
 
         msg = 'Vertex indices reference non-existing coordinate sets'
-        assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
+        assert max(max(self.triangles)) <= self.nodes.shape[0], msg
 
 
         # FIXME: Maybe move to statistics?
         # Or use with get_extent
-        xy_extent = [ min(self.coordinates[:,0]), min(self.coordinates[:,1]) ,
-                      max(self.coordinates[:,0]), max(self.coordinates[:,1]) ]
+        xy_extent = [ min(self.nodes[:,0]), min(self.nodes[:,1]) ,
+                      max(self.nodes[:,0]), max(self.nodes[:,1]) ]
 
         self.xy_extent = array(xy_extent, Float)
@@ -152 +165 @@
             if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
            
-
-            x0 = V[i, 0]; y0 = V[i, 1]
-            x1 = V[i, 2]; y1 = V[i, 3]
-            x2 = V[i, 4]; y2 = V[i, 5]
+            x0, y0 = V[3*i, :]
+            x1, y1 = V[3*i+1, :]
+            x2, y2 = V[3*i+2, :]            
 
             # Area
@@ -210 +222 @@
     def __repr__(self):
         return 'Mesh: %d vertices, %d triangles'\
-               %(self.coordinates.shape[0], len(self))
+               %(self.nodes.shape[0], len(self))
 
     def get_normals(self):
@@ -227 +239 @@
 
 
-
-    def get_vertex_coordinates(self, unique=False, obj=False, absolute=False):
-        """Return all vertex coordinates.
-        Return all vertex coordinates for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
-
-        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
-        FIXME, we might make that the default.
-        FIXME Maybe use keyword: continuous = False for this condition?
-        FIXME - Maybe use something referring to unique vertices?
+    def get_nodes(self, absolute=False):
+        """Return all node coordinates ordered in an Nx2 array.
+        This is the same format they were provided in the constructor
+        i.e. without any duplication.
 
         Boolean keyword argument absolute determines whether coordinates
         are to be made absolute by taking georeference into account
         Default is False as many parts of ANUGA expects relative coordinates.
-        (To see which, switch to default absolute=True and run tests).
+        (To see which, switch to default absolute=True and run tests).        
         """
 
         N = self.number_of_full_nodes
-
-        if unique is True:
-            V = self.coordinates[:N,:]
-            if absolute is True:
-                if not self.geo_reference.is_absolute():
-                    V = self.geo_reference.get_absolute(V)
-
-            return V
-
+        V = self.nodes[:N,:]
+        if absolute is True:
+            if not self.geo_reference.is_absolute():
+                V = self.geo_reference.get_absolute(V)
                 
-
+        return V
+        
+        
+
+    def get_vertex_coordinates(self, unique=False, absolute=False):
+        """Return all vertex coordinates.
+        Return all vertex coordinates for all triangles as a 3*N x 2 array
+        where the jth vertex of the ith triangle is located in row 3*i+j.
+
+        Boolean keyword unique will cause the points to be returned as
+        they were provided in the constructor i.e. without any duplication
+        in an N x 2 array.
+
+        """
+
+        if unique is True:        
+            return self.get_nodes(absolute)
+
+            
         V = self.vertex_coordinates
         if absolute is True:
             if not self.geo_reference.is_absolute():
+                V = self.geo_reference.get_absolute(V)
             
-                V0 = self.geo_reference.get_absolute(V[:N,0:2])
-                V1 = self.geo_reference.get_absolute(V[:N,2:4])
-                V2 = self.geo_reference.get_absolute(V[:N,4:6])
-
-                # This does double the memory need 
-                V = concatenate( (V0, V1, V2), axis=1 )
-
-                
-        if obj is True:
-            N = V.shape[0]
-            return reshape(V, (3*N, 2))
-        else:    
-            return V
+        return V
+
 
 
@@ -281 +290 @@
 
         V = self.get_vertex_coordinates(absolute=absolute)
-        return V[i, 2*j:2*j+2]
-    
-        ##return self.vertex_coordinates[i, 2*j:2*j+2]
+        return V[3*i+j, :]
 
 
     def compute_vertex_coordinates(self):
-        """Return vertex coordinates for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
-        """
-
-        #FIXME (Ole): Perhaps they should be ordered as in obj files??
-        #See quantity.get_vertex_values
-        #FIXME (Ole) - oh yes they should
+        """Return all vertex coordinates for all triangles as a 3*N x 2 array
+        where the jth vertex of the ith triangle is located in row 3*i+j.
+
+        This function is used to precompute this important structure. Use
+        get_vertex coordinates to retrieve the points.
+        """
 
         N = self.number_of_triangles
-        vertex_coordinates = zeros((N, 6), Float)
+        vertex_coordinates = zeros((3*N, 2), Float)
 
         for i in range(N):
             for j in range(3):
                 k = self.triangles[i,j]  #Index of vertex 0
-                v_k = self.coordinates[k]
-                vertex_coordinates[i, 2*j+0] = v_k[0]
-                vertex_coordinates[i, 2*j+1] = v_k[1]
+                v_k = self.nodes[k]
+
+                vertex_coordinates[3*i+j,:] = v_k
+
 
         return vertex_coordinates
+
 
     def get_vertices(self, indices=None):
@@ -358 +366 @@
 
         Preconditions:
-          self.coordinates and self.triangles are defined
+          self.nodes and self.triangles are defined
 
         Postcondition:
@@ -364 +372 @@
         """
 
-        vertexlist = [None]*len(self.coordinates)
+        vertexlist = [None]*len(self.nodes)
         for i in range(self.number_of_triangles):
 

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -220 +220 @@
         for i, (vol_id, edge_id) in enumerate(boundary_keys):
 
-            x0 = V[vol_id, 0]; y0 = V[vol_id, 1]
-            x1 = V[vol_id, 2]; y1 = V[vol_id, 3]
-            x2 = V[vol_id, 4]; y2 = V[vol_id, 5]
-
+            base_index = 3*vol_id
+            x0, y0 = V[base_index, :]
+            x1, y1 = V[base_index+1, :]
+            x2, y2 = V[base_index+2, :]
+            
             #Compute midpoints
             if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -113 +113 @@
             if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
 
-            x0 = V[i, 0]; y0 = V[i, 1]
-            x1 = V[i, 2]; y1 = V[i, 3]
-            x2 = V[i, 4]; y2 = V[i, 5]
+            x0, y0 = V[3*i, :]
+            x1, y1 = V[3*i+1, :]
+            x2, y2 = V[3*i+2, :]                        
+
+            #x0 = V[i, 0]; y0 = V[i, 1]
+            #x1 = V[i, 2]; y1 = V[i, 3]
+            #x2 = V[i, 4]; y2 = V[i, 5]
 
             #Compute centroid
@@ -596 +600 @@
         #Check each triangle
         for i in range(N):
+
+            x0, y0 = V[3*i, :]
+            x1, y1 = V[3*i+1, :]
+            x2, y2 = V[3*i+2, :]
             
-            x0 = V[i, 0]; y0 = V[i, 1]
-            x1 = V[i, 2]; y1 = V[i, 3]
-            x2 = V[i, 4]; y2 = V[i, 5]
-
             #Check that area hasn't been compromised
             area = self.areas[i]

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -16 +16 @@
 
 from Numeric import array, zeros, Float, less, concatenate, NewAxis, argmax, allclose
-from anuga.utilities.numerical_tools import ensure_numeric
+from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
 
 class Quantity:
@@ -558 +558 @@
             is_subset = True
 
+
+        # FIXME (Ole): Now we can compute the arrays once and for all
+        # for both centroids and vertices and the use set_values_from_array
+
         if location == 'centroids':
-            P = take(self.domain.centroid_coordinates, indices)
+            V = take(self.domain.centroid_coordinates, indices)
             if is_subset:
-                self.set_values(f(P[:,0], P[:,1]),
+                self.set_values(f(V[:,0], V[:,1]),
                                 location = location,
                                 indices = indices)
             else:
-                self.set_values(f(P[:,0], P[:,1]), location = location)
+                self.set_values(f(V[:,0], V[:,1]), location = location)
         elif location == 'vertices':
-            P = self.domain.vertex_coordinates
+            V = self.domain.get_vertex_coordinates()
+
+            x = V[:,0]; y = V[:,1];                     
+            values = f(x, y)
+
+            if is_scalar(values):
+                # Function returned a constant value
+                self.set_values_from_constant(values,
+                                              location, indices, verbose)
+                return
+                
+            
             if is_subset:
                 #Brute force
-                for e in indices:
-                    for i in range(3):
-                        self.vertex_values[e,i] = f(P[e,2*i], P[e,2*i+1])
+                for i in indices:
+                    for j in range(3):
+                        self.vertex_values[i,j] = values[3*i+j]
             else:
-                for i in range(3):
-                    self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1])
+                for j in range(3):
+                    self.vertex_values[:,j] = values[j::3] 
         else:
             raise 'Not implemented: %s' %location
@@ -625 +640 @@
             raise ms
 
-        coordinates = self.domain.coordinates
-        triangles = self.domain.triangles
+        coordinates = self.domain.get_nodes()
+        triangles = self.domain.triangles      #FIXME
 
 
@@ -910 +925 @@
         elif location == 'unique vertices':
             if (indices ==  None):
-                indices=range(self.domain.coordinates.shape[0])
+                indices=range(self.domain.number_of_nodes)
             vert_values = []
             #Go through list of unique vertices
@@ -958 +973 @@
 
         if indices == None:
-            assert A.shape[0] == self.domain.coordinates.shape[0]
+            assert A.shape[0] == self.domain.get_nodes().shape[0]
             vertex_list = range(A.shape[0])
         else:
@@ -1085 +1100 @@
 
             if xy is True:
-                X = self.domain.coordinates[:,0].astype(precision)
-                Y = self.domain.coordinates[:,1].astype(precision)
+                X = self.domain.get_nodes()[:,0].astype(precision)
+                Y = self.domain.get_nodes()[:,1].astype(precision)
 
                 return X, Y, A, V
@@ -1389 +1404 @@
     from Numeric import zeros, Float
 
-    N = len(quantity.domain)
+    N = quantity.domain.number_of_nodes
 
     beta_w = quantity.domain.beta_w

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -29 +29 @@
 
 
-        assert allclose(domain.get_vertex_coordinates(unique=True), domain.coordinates)
+        assert allclose(domain.get_vertex_coordinates(unique=True),
+                        domain.nodes)
 
         #assert allclose(domain.get_vertex_coordinates(), ...TODO
@@ -51 +52 @@
         assert domain.get_vertices([0,4]) == [domain.triangles[0],
                                               domain.triangles[4]]
+        
     def test_areas(self):
         from mesh_factory import rectangular

anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -76 +76 @@
 
         #Vertex coordinates
+        #V = mesh.get_vertex_coordinates()
+        #assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0])
+        
+
         V = mesh.get_vertex_coordinates()
-        assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0])
-
-        V = mesh.get_vertex_coordinates(obj=True)
         assert allclose(V, [ [0.0, 0.0],
                              [4.0, 0.0],
@@ -104 +105 @@
 
         V = mesh.get_vertex_coordinates()
-
-        x0 = V[0,0]
-        y0 = V[0,1]
-        x1 = V[0,2]
-        y1 = V[0,3]
-        x2 = V[0,4]
-        y2 = V[0,5]
+        x0 = V[0, 0]; y0 = V[0, 1]
+        x1 = V[1, 0]; y1 = V[1, 1]
+        x2 = V[2, 0]; y2 = V[2, 1]
+        #x0 = V[0,0]
+        #y0 = V[0,1]
+        #x1 = V[0,2]
+        #y1 = V[0,3]
+        #x2 = V[0,4]
+        #y2 = V[0,5]
 
         m0 = [(x1 + x2)/2, (y1 + y2)/2]
@@ -172 +175 @@
 
         V = mesh.get_vertex_coordinates()
-
-        x0 = V[0,0]
-        y0 = V[0,1]
-        x1 = V[0,2]
-        y1 = V[0,3]
-        x2 = V[0,4]
-        y2 = V[0,5]
+        x0 = V[0, 0]; y0 = V[0, 1]
+        x1 = V[1, 0]; y1 = V[1, 1]
+        x2 = V[2, 0]; y2 = V[2, 1]        
+
+        #x0 = V[0,0]
+        #y0 = V[0,1]
+        #x1 = V[0,2]
+        #y1 = V[0,3]
+        #x2 = V[0,4]
+        #y2 = V[0,5]
 
         m0 = [(x1 + x2)/2, (y1 + y2)/2]

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -395 +395 @@
 
 
-        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
+        answer = linear_function(quantity.domain.get_vertex_coordinates())
         #print quantity.vertex_values, answer
         assert allclose(quantity.vertex_values.flat, answer)
@@ -401 +401 @@
 
         #Now try by setting the same values directly
-        vertex_attributes = fit_to_mesh(quantity.domain.coordinates,
-                                        quantity.domain.triangles,
+        vertex_attributes = fit_to_mesh(quantity.domain.get_nodes(),
+                                        quantity.domain.triangles, #FIXME
                                         data_points,
                                         z,
@@ -514 +514 @@
         quantity.set_values(filename = ptsfile,
                             attribute_name = att, alpha = 0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
+        answer = linear_function(quantity.domain.get_vertex_coordinates())
 
         #print quantity.vertex_values.flat
@@ -592 +592 @@
         quantity.set_values(filename = ptsfile,
                             attribute_name = att, alpha = 0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True) - [x0, y0])
+        answer = linear_function(quantity.domain.get_vertex_coordinates() - [x0, y0])
 
         assert allclose(quantity.vertex_values.flat, answer)
@@ -666 +666 @@
         quantity.set_values(filename = ptsfile,
                             attribute_name = att, alpha = 0)
-        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True) + [x0, y0])
+        answer = linear_function(quantity.domain.get_vertex_coordinates() + [x0, y0])
 
 

anuga_core/source/anuga/fit_interpolate/fit.py

@@ -93 +93 @@
                  max_vertices_per_cell)
         
-        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (vertices)
+        m = self.mesh.number_of_nodes # Nbr of basis functions (vertices)
         
         self.AtA = None
@@ -152 +152 @@
         #"element stiffness matrices:
 
-        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
+        m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
 
         self.D = Sparse(m,m)
@@ -233 +233 @@
         if self.AtA == None:
             # AtA and Atz need ot be initialised.
-            m = self.mesh.coordinates.shape[0] #Nbr of vertices
+            m = self.mesh.number_of_nodes
             if len(z.shape) > 1:
                 att_num = z.shape[1]
@@ -316 +316 @@
 
         #Check sanity
-        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
+        m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
         n = self.point_count
         if n<m and self.alpha == 0.0:

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -257 +257 @@
             print '\n WARNING: No points within the mesh! \n'
             
-        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
-        n = point_coordinates.shape[0]     #Nbr of data points
+        m = self.mesh.number_of_nodes  # Nbr of basis functions (1/vertex)
+        n = point_coordinates.shape[0] # Nbr of data points
 
         if verbose: print 'Number of datapoints: %d' %n

anuga_core/source/anuga/utilities/numerical_tools.py

@@ -61 +61 @@
     if x < 0: return -1
     if x == 0: return 0    
-    
+
+
+def is_scalar(x):
+    """True if x is a scalar (constant numeric value)
+    """
+
+    from types import IntType, FloatType
+    if type(x) in [IntType, FloatType]:
+        return True
+    else:
+        return False
 
 def angle(v1, v2=None):

anuga_core/source/anuga/utilities/quad.py

@@ -139 +139 @@
         if len(args) == 2:
             point_id = int(args[1])
-            x = self.__class__.mesh.coordinates[point_id][0]
-            y = self.__class__.mesh.coordinates[point_id][1]
+            x, y = self.__class__.mesh.get_nodes()[point_id]
 
             #print point_id, x, y
@@ -394 +393 @@
     #Make root cell
     #print mesh.coordinates
-    
-    xmin = min(mesh.coordinates[:,0])
-    xmax = max(mesh.coordinates[:,0])
-    ymin = min(mesh.coordinates[:,1])
-    ymax = max(mesh.coordinates[:,1])
+
+    nodes = mesh.get_nodes()
+    xmin = min(nodes[:,0])
+    xmax = max(nodes[:,0])
+    ymin = min(nodes[:,1])
+    ymax = max(nodes[:,1])
 
     
@@ -425 +425 @@
     
     #Insert indices of all vertices
-    root.insert( range(len(mesh.coordinates)) )
+    root.insert( range(mesh.number_of_nodes) )
 
     #Build quad tree and return

anuga_core/source/anuga/utilities/test_quad.py

@@ -120 +120 @@
 
         Q = build_quadtree(self.mesh)
+        #Q.show()
+        #print Q.count()
         assert Q.count() == 8
 
-        #Q.show()
+
 
         result = Q.search(3, 105)