Changeset 7710


Ignore:
Timestamp:
May 7, 2010, 12:25:58 PM (15 years ago)
Author:
James Hudson
Message:

Refactored quadtree interface.

Location:
anuga_core/source/anuga
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

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

    r7707 r7710  
    5757                 mesh_origin=None,
    5858                 alpha = None,
    59                  verbose=False,
    60                  max_vertices_per_cell=None):
     59                 verbose=False):
    6160
    6261
    … …  
    7978              If specified vertex coordinates are assumed to be
    8079              relative to their respective origins.
    81 
    82           max_vertices_per_cell: Number of vertices in a quad tree cell
    83           at which the cell is split into 4.
    8480
    8581          Note: Don't supply a vertex coords as a geospatial object and
    … …  
    10399                 mesh,
    104100                 mesh_origin,
    105                  verbose,
    106                  max_vertices_per_cell)
     101                 verbose)
    107102       
    108103        m = self.mesh.number_of_nodes # Nbr of basis functions (vertices)

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

    r7690 r7710  
    5454                 mesh=None,
    5555                 mesh_origin=None,
    56                  verbose=False,
    57                  max_vertices_per_cell=None):
     56                 verbose=False):
    5857
    5958
    … …  
    8281              relative to their respective origins.
    8382
    84           max_vertices_per_cell: Number of vertices in a quad tree cell
    85           at which the cell is split into 4.
    86 
    8783          Note: Don't supply a vertex coords as a geospatial object and
    8884              a mesh origin, since geospatial has its own mesh origin.
    8985        """
    9086        global build_quadtree_time
    91         if max_vertices_per_cell == None:
    92             max_vertices_per_cell = MAX_VERTICES_PER_CELL
    9387        if mesh is None:
    9488            if vertex_coordinates is not None and  triangles is not None:
    … …  
    113107            #This stores indices of vertices
    114108            t0 = time.time()
    115             self.root = build_quadtree(self.mesh,
    116                                        max_points_per_cell = max_vertices_per_cell)
     109            self.root = build_quadtree(self.mesh)
    117110       
    118111            build_quadtree_time =  time.time()-t0

  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r7707 r7710  
    5757# @param interpolation_points Points to interpolate to.
    5858# @param mesh_origin A geo_ref object or 3-tuples of UTMzone, easting, northing.
    59 # @param max_vertices_per_cell Max number of vertices before splitting cell.
    6059# @param start_blocking_len Block if # of points greater than this.
    6160# @param use_cache If True, cache.
    … …  
    6665                interpolation_points,
    6766                mesh_origin=None,
    68                 max_vertices_per_cell=None,
    6967                start_blocking_len=500000,
    7068                use_cache=False,
    … …  
    103101                 relative to their respective origins.
    104102
    105     max_vertices_per_cell: Number of vertices in a quad tree cell
    106                            at which the cell is split into 4.
    107 
    108103                           Note: Don't supply a vertex coords as a geospatial
    109104                           object and a mesh origin, since geospatial has its
    … …  
    136131            ensure_numeric(triangles))
    137132    kwargs = {'mesh_origin': mesh_origin,
    138               'max_vertices_per_cell': max_vertices_per_cell,
    139133              'verbose': verbose}
    140134
    … …  
    175169                 triangles,
    176170                 mesh_origin=None,
    177                  verbose=False,
    178                  max_vertices_per_cell=None):
     171                 verbose=False):
    179172
    180173        """ Build interpolation matrix mapping from
    … …  
    209202                                triangles=triangles,
    210203                                mesh_origin=mesh_origin,
    211                                 verbose=verbose,
    212                                 max_vertices_per_cell=max_vertices_per_cell)
     204                                verbose=verbose)
    213205
    214206        # Initialise variables

  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r7685 r7710  
    307307
    308308        data = [ [4,4] ]
    309         interp = Interpolate(points, triangles,
    310                                max_vertices_per_cell = 4)
     309        interp = Interpolate(points, triangles)
    311310        #print "PDSG - interp.get_A()", interp.get_A()
    312311        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,

  • anuga_core/source/anuga/fit_interpolate/test_search_functions.py

    r7707 r7710  
    5050        mesh.check_integrity()
    5151
    52         root = build_quadtree(mesh, max_points_per_cell = 1)
     52        root = build_quadtree(mesh)
    5353        set_last_triangle()
    5454
    … …  
    5858        found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, [0, 0])
    5959        assert found is True
    60                
    61                
     60       
     61       
    6262    def test_small(self):
    6363        """test_small: Two triangles
    … …  
    7070        mesh.check_integrity()
    7171
    72         root = build_quadtree(mesh, max_points_per_cell = 1)
     72        root = build_quadtree(mesh)
    7373        set_last_triangle()
    7474
    … …  
    7676        found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
    7777        assert k == 1 # Triangle one
    78         assert found is True           
    79                
     78        assert found is True       
     79       
    8080    def test_bigger(self):
    8181        """test_bigger
    … …  
    9090        mesh.check_integrity()
    9191
    92         root = build_quadtree(mesh, max_points_per_cell = 4)
     92        root = build_quadtree(mesh)
    9393        set_last_triangle()
    9494
    … …  
    9898           
    9999            found, s0, s1, s2, k = search_tree_of_vertices(root, mesh,
    100                                                            ensure_numeric(x))                                                              
    101                                                                                                                    
     100                                                           ensure_numeric(x))                                   
     101                                                           
    102102            if k >= 0:
    103103                V = mesh.get_vertex_coordinates(k) # nodes for triangle k
    … …  
    121121
    122122       
    123         for m in range(8):
    124             root = build_quadtree(mesh, max_points_per_cell = m)
    125             set_last_triangle()
    126             #print m, root.show()
    127 
    128             for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7],
    129                       [0.1,0.9], [0.4,0.6], [0.9,0.1],
    130                       [10, 3]]:
    131                
    132                 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
    133 
    134                 if k >= 0:
    135                     V = mesh.get_vertex_coordinates(k) # nodes for triangle k
    136                     assert is_inside_triangle(x, V, closed=True)
    137                     assert is_inside_polygon(x, V)
    138                     assert found is True
    139                 else:
    140                     assert found is False               
    141 
     123
     124        root = build_quadtree(mesh)
     125        set_last_triangle()
     126        #print m, root.show()
     127
     128        for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7],
     129                  [0.1,0.9], [0.4,0.6], [0.9,0.1],
     130                  [10, 3]]:
    142131           
    143                 if m == k == 0: return   
     132            found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
     133
     134            if k >= 0:
     135                V = mesh.get_vertex_coordinates(k) # nodes for triangle k
     136                assert is_inside_triangle(x, V, closed=True)
     137                assert is_inside_polygon(x, V)
     138                assert found is True
     139            else:
     140                assert found is False               
     141
     142       
     143            if k == 0: return   
    144144
    145145    def test_underlying_function(self):
    … …  
    150150        mesh = Mesh(points, vertices, boundary)
    151151
    152         root = build_quadtree(mesh, max_points_per_cell = 4)
     152        root = build_quadtree(mesh)
    153153        set_last_triangle()
    154154
    … …  
    157157
    158158        triangles = _trilist_from_indices(mesh, root.search(x[0], x[1]))
    159        
     159   
    160160        found, sigma0, sigma1, sigma2, k = \
    161161               _search_triangles_of_vertices(triangles, x)

  • anuga_core/source/anuga/utilities/quad.py

    r7709 r7710  
    99import string, types, sys
    1010import anuga.utilities.log as log
     11from anuga.utilities.aabb import AABB
    1112
    12 # Allow children to be slightly bigger than their parents to prevent straddling of a boundary
    13 SPLIT_BORDER_RATIO    = 0.55
    14 
    15 class AABB:
    16     """Axially-aligned bounding box class.
    17     """
    18    
    19     def __init__(self, xmin, xmax, ymin, ymax):
    20         self.xmin = xmin   
    21         self.xmax = xmax
    22         self.ymin = ymin   
    23         self.ymax = ymax
    24 
    25     def __repr__(self):
    26         return '(xmin:%f, xmax:%f, ymin:%f, ymax:%f)' \
    27                % (round(self.xmin,1), round(self.xmax,1), round(self.ymin,1), round(self.ymax, 1))
    28        
    29     def size(self):
    30         """return size as (w,h)"""
    31         return self.xmax - self.xmin, self.ymax - self.ymin
    32        
    33     def split(self, border=SPLIT_BORDER_RATIO):
    34         """Split along shorter axis.
    35            return 2 subdivided AABBs.
    36         """
    37        
    38         width, height = self.size()
    39         assert width >= 0 and height >= 0
    40        
    41         if (width > height):
    42             # split vertically
    43             return AABB(self.xmin, self.xmin+width*border, self.ymin, self.ymax), \
    44                    AABB(self.xmax-width*border, self.xmax, self.ymin, self.ymax)
    45         else:
    46             # split horizontally       
    47             return AABB(self.xmin, self.xmax, self.ymin, self.ymin+height*border), \
    48                    AABB(self.xmin, self.xmax, self.ymax-height*border, self.ymax)   
    49    
    50     def is_trivial_in(self, test):
    51         if (test.xmin < self.xmin) or (test.xmax > self.xmax):
    52             return False       
    53         if (test.ymin < self.ymin) or (test.ymax > self.ymax):
    54             return False       
    55         return True
    56  
    57     def contains(self, x, y):
    58         return (self.xmin <= x <= self.xmax) and (self.ymin <= y <= self.ymax)
    5913           
    6014class Cell(TreeNode):
    … …  
    6317    One cell in the plane delimited by southern, northern,
    6418    western, eastern boundaries.
    65 
    66     Public Methods:
    67         insert(point)
    68         search(x, y)
    69         split()
    70         store()
    71         retrieve()
    72         count()
    7319    """
    7420 
    … …  
    8834    def __repr__(self):
    8935        str = '%s: leaves: %d' \
    90                % (self.name , len(self.leaves)) 
     36               % (self.name , len(self.leaves))   
    9137        if self.children:
    9238            str += ', children: %d' % (len(self.children))
    … …  
    211157#from anuga.pmesh.mesh import Mesh
    212158   
    213 def build_quadtree(mesh, max_points_per_cell = 4):
     159def build_quadtree(mesh):
    214160    """Build quad tree for mesh.
    215161
    216     All vertices in mesh are stored in quadtree and a reference
     162    All vertex indices in mesh are stored in a quadtree and a reference
    217163    to the root is returned.
    218164    """
    219 
    220 
    221     #Make root cell
    222     #print mesh.coordinates
    223 
    224     xmin, xmax, ymin, ymax = mesh.get_extent(absolute=True)
    225    
    226     # Ensure boundary points are fully contained in region
    227     # It is a property of the cell structure that
    228     # points on xmax or ymax of any given cell
    229     # belong to the neighbouring cell.
    230     # Hence, the root cell needs to be expanded slightly
    231     ymax += (ymax-ymin)/10
    232     xmax += (xmax-xmin)/10
    233 
    234     # To avoid round off error
    235     ymin -= (ymax-ymin)/10
    236     xmin -= (xmax-xmin)/10   
    237 
    238     #print "xmin", xmin
    239     #print "xmax", xmax
    240     #print "ymin", ymin
    241     #print "ymax", ymax
    242    
    243     root = Cell(AABB(xmin, xmax, ymin, ymax))
     165    extents = AABB(*mesh.get_extent(absolute=True))   
     166    extents.grow(1.001) # To avoid round off error
     167    root = Cell(extents)
    244168   
    245169    N = len(mesh)
    … …  
    247171    # Get x,y coordinates for all vertices for all triangles
    248172    V = mesh.get_vertex_coordinates(absolute=True)
    249        
     173   
    250174    # Check each triangle
    251175    for i in range(N):
Note: See TracChangeset for help on using the changeset viewer.