Changeset 4655


Ignore:
Timestamp:
Aug 6, 2007, 1:53:26 PM (17 years ago)
Author:
duncan
Message:

When fitting, cell data - triangle vertices and norms - are calculated the first time a point is looked for in a cell. This is to speed things up. Also, the old point info is deleted to save memory.

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/search_functions.py

    r4654 r4655  
    3737    k = -10.0
    3838           
    39     #Find vertices near x
    40     candidate_vertices = root.search(x[0], x[1])
     39    #Get triangles in the cell that the point is in.
     40    # Triangle is a list, first element triangle_id,
     41    # second element the triangle
     42    triangles = root.search(x[0], x[1])
    4143    is_more_elements = True
    4244
    4345    element_found, sigma0, sigma1, sigma2, k = \
    4446                   _search_triangles_of_vertices(mesh,
    45                                                  candidate_vertices, x)
     47                                                 triangles, x)
    4648    while not element_found and is_more_elements:
    47         candidate_vertices, branch = root.expand_search()
     49        triangles, branch = root.expand_search()
    4850        if branch == []:
    4951            # Searching all the verts from the root cell that haven't
    5052            # been searched.  This is the last try
    5153            element_found, sigma0, sigma1, sigma2, k = \
    52                            _search_triangles_of_vertices(mesh,
    53                                                          candidate_vertices, x)
     54                           _search_triangles_of_vertices(mesh,triangles, x)
    5455            is_more_elements = False
    5556        else:
    5657            element_found, sigma0, sigma1, sigma2, k = \
    57                        _search_triangles_of_vertices(mesh,
    58                                                      candidate_vertices, x)
     58                       _search_triangles_of_vertices(mesh,triangles, x)
    5959
    6060    return element_found, sigma0, sigma1, sigma2, k
    6161
    6262
    63 def _search_triangles_of_vertices(mesh, candidate_vertices, x):
     63def _search_triangles_of_vertices(mesh, triangles, x):
    6464    """Search for triangle containing x amongs candidate_vertices in mesh
    6565
     
    8282   
    8383    #For all vertices in same cell as point x
    84     for v in candidate_vertices:
    85        
    86         #FIXME (DSG-DSG): this catches verts with no triangle.
    87         #Currently pmesh is producing these.
    88         #this should be stopped,
    89 
    90         if mesh.number_of_triangles_per_node[v] == 0:
    91             continue
    92        
    93         # Get all triangles which has v as a vertex
    94         # The list has elements (triangle, vertex), but only the
    95         # first component will be used here
    96         triangle_list = mesh.get_triangles_and_vertices_per_node(node=v)
    97 
     84    for k, tri_verts_norms in triangles:
     85        tri = tri_verts_norms[0]
     86        n0, n1, n2 = tri_verts_norms[1]
     87        # k is the triangle index
     88        # tri is a list of verts (x, y), representing a tringle
    9889        # Find triangle that contains x (if any) and interpolate
    99        
    100         for k, _ in triangle_list:
    101             element_found, sigma0, sigma1, sigma2, k =\
    102                            find_triangle_compute_interpolation(mesh,
    103                                                                k,
    104                                                                x)
    105 
    106             if element_found is True:
    107                 # Don't look for any other triangles in the triangle list
    108                 break
    109 
     90        element_found, sigma0, sigma1, sigma2 =\
     91                       find_triangle_compute_interpolation(tri, n0, n1, n2, x)
    11092        if element_found is True:
    111             # Don't look for any other triangle_lists from the
    112             # candidate_vertices
     93            # Don't look for any other triangles in the triangle list
    11394            break
    114        
    11595    return element_found, sigma0, sigma1, sigma2, k
    11696
    11797
    11898           
    119 def find_triangle_compute_interpolation(mesh, k, x):
     99def find_triangle_compute_interpolation(triangle, n0, n1, n2, x):
    120100    """Compute linear interpolation of point x and triangle k in mesh.
    121101    It is assumed that x belongs to triangle k.
     
    123103
    124104    # Get the three vertex_points of candidate triangle k
    125     xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)
     105    xi0, xi1, xi2 = triangle
    126106
    127107    # this is where we can call some fast c code.
     
    138118   
    139119    if  x[0] > xmax + epsilon:
    140         return False,0,0,0,0
     120        return False,0,0,0
    141121    if  x[0] < xmin - epsilon:
    142         return False,0,0,0,0
     122        return False,0,0,0
    143123    if  x[1] > ymax + epsilon:
    144         return False,0,0,0,0
     124        return False,0,0,0
    145125    if  x[1] < ymin - epsilon:
    146         return False,0,0,0,0
     126        return False,0,0,0
    147127   
    148128    # Get the three normals
    149     n0 = mesh.get_normal(k, 0)
    150     n1 = mesh.get_normal(k, 1)
    151     n2 = mesh.get_normal(k, 2)           
    152        
     129    #n0 = norms[0] 
     130    #n1 = norms[1]
     131    #n2 = norms[2]
    153132       
    154133    # Compute interpolation
     
    171150    else:
    172151        element_found = False
    173     return element_found, sigma0, sigma1, sigma2, k
     152    return element_found, sigma0, sigma1, sigma2
  • anuga_core/source/anuga/fit_interpolate/test_search_functions.py

    r4648 r4655  
    1212from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
    1313from anuga.utilities.polygon import is_inside_polygon
    14 from anuga.utilities.quad import build_quadtree
    15 
    16 
    17 
     14from anuga.utilities.quad import build_quadtree, Cell
    1815
    1916
     
    161158                assert found is True
    162159            else:
    163                 assert found is False   
     160                assert found is False
     161
     162               
     163
     164    def expanding_search(self):
     165        """test_larger mesh and different quad trees
     166        """
     167       
     168        p0 = [2,1]
     169        p1 = [4,1]
     170        p2 = [4.,4]
     171        p3 = [2,4]
     172        p4 = [5,4]
     173
     174        p5 = [-1,-1]
     175        p6 = [1,-1]
     176        p7 = [1,1]
     177        p8 = [-1,1]
     178
     179        points = [p0,p1,p2, p3,p4,p5,p6,p7,p8]
     180        #
     181        vertices = [[0,1,2],[0,2,3],[1,4,2],[5,6,7], [5,7,8]]
     182        mesh = Mesh(points, vertices)
     183
     184        # Don't do this, want to control the mx and mins
     185        #root = build_quadtree(mesh, max_points_per_cell=4)
     186   
     187        #Initialise
     188        Cell.initialise(mesh)
     189
     190        root = Cell(-3, 9, -3, 9,
     191                    max_points_per_cell = 4)
     192        #Insert indices of all vertices
     193        root.insert( range(mesh.number_of_nodes) )
     194
     195        #Build quad tree and return
     196        root.split()
     197       
     198        # One point
     199        #x = [3.5, 1.5]
     200        x = [2.5, 1.5]
     201        element_found, sigma0, sigma1, sigma2, k = \
     202                       search_tree_of_vertices(root, mesh, x)
     203        # One point
     204        x = [3.00005, 2.999994]
     205        element_found, sigma0, sigma1, sigma2, k = \
     206                       search_tree_of_vertices(root, mesh, x)
     207        assert element_found is True
     208        assert k == 1
     209       
     210
    164211#-------------------------------------------------------------
    165212if __name__ == "__main__":
    166213    suite = unittest.makeSuite(Test_search_functions,'test')
     214    #suite = unittest.makeSuite(Test_search_functions,'expanding_search')
    167215    runner = unittest.TextTestRunner(verbosity=1)
    168216    runner.run(suite)
  • anuga_core/source/anuga/fit_interpolate/ticket178_benchmark.py

    r4654 r4655  
    1414ofile = 'lbm_resultsII.csv'
    1515delimiter = ','
    16 run_profile = True #False #True
     16run_profile = False #True
    1717is_fit_list = [True, False]
    1818num_of_points_list = [3, 200, 600, 2000, 6000, 10000, 20000]
  • anuga_core/source/anuga/utilities/quad.py

    r4654 r4655  
    7272       
    7373 
    74     def search(self, x, y):
    75     #def search_new(self, x, y):
     74    def search(self, x, y, get_vertices=False):
    7675        """Find all point indices sharing the same cell as point (x, y)
    7776        """
     
    8483                    brothers.remove(child)
    8584                    branch.append(brothers)
    86                     points, branch = child.search_branch(x,y, branch)
     85                    points, branch = child.search_branch(x,y, branch,
     86                                                  get_vertices=get_vertices)
    8787        else:
    8888            # Leaf node: Get actual waypoints
    89             points = self.retrieve()
    90 
     89            points = self.retrieve(get_vertices=get_vertices)
    9190        self.branch = branch   
    9291        return points
    9392
    9493
    95     def search_branch(self, x, y, branch):
     94    def search_branch(self, x, y, branch, get_vertices=False):
    9695        """Find all point indices sharing the same cell as point (x, y)
    9796        """
     
    103102                    brothers.remove(child)
    104103                    branch.append(brothers)
    105                     points, branch = child.search_branch(x,y, branch)
     104                    points, branch = child.search_branch(x,y, branch,
     105                                                  get_vertices=get_vertices)
    106106                   
    107107        else:
    108108            # Leaf node: Get actual waypoints
    109             points = self.retrieve()     
     109            points = self.retrieve(get_vertices=get_vertices)     
    110110        return points, branch
    111111
    112112
    113     def expand_search(self):
     113    def expand_search(self, get_vertices=False):
    114114        """Find all point indices 'up' one cell from the last search
    115115        """
     116       
    116117        points = []
    117118        if self.branch == []:
     
    121122            for cell in three_cells:
    122123                #print "cell ", cell.show()
    123                 points += cell.retrieve()
     124                points += cell.retrieve(get_vertices=get_vertices)
    124125        return points, self.branch
    125126
     
    200201
    201202
    202     def retrieve(self):
     203    def retrieve_triangles(self):
     204        """return a list of lists. For the inner lists,
     205        The first element is the triangle index,
     206        the second element is a list.for this list
     207           the first element is a list of three (x, y) vertices,
     208           the following elements are the three triangle normals.
     209
     210        This info is used in searching for a triangle that a point is in.
     211
     212        Post condition
     213        No more points can be added to the quad tree, since the
     214        points data structure is removed.
     215        """
     216        # FIXME Tidy up the structure that is returned.
     217        # if the triangles att has been made
     218        # return it.
     219        if not hasattr(self,'triangles'):
     220            # use a dictionary to remove duplicates
     221            triangles = {}
     222            verts = self.retrieve_vertices()
     223            # print "verts", verts
     224            for vert in verts:
     225                triangle_list = self.__class__.mesh.get_triangles_and_vertices_per_node(vert)
     226                for k, _ in triangle_list:
     227                    if not triangles.has_key(k):
     228                        # print 'k',k
     229                        tri = self.__class__.mesh.get_vertex_coordinates(k)
     230                        n0 = self.__class__.mesh.get_normal(k, 0)
     231                        n1 = self.__class__.mesh.get_normal(k, 1)
     232                        n2 = self.__class__.mesh.get_normal(k, 2)
     233                        triangles[k]=(tri, (n0, n1, n2))
     234            self.triangles = triangles.items()
     235            # Delete the old cell data structure to save memory
     236            del self.points
     237        return self.triangles
     238           
     239    def retrieve_vertices(self):
     240         return self.points
     241
     242
     243    def retrieve(self, get_vertices=True):
    203244         objects = []
    204245         if self.children is None:
    205              objects = self.points
     246             if get_vertices is True:
     247                 objects = self.retrieve_vertices()
     248             else:
     249                 objects =  self.retrieve_triangles()
    206250         else: 
    207251             for child in self:
    208                  objects += child.retrieve()
    209          return objects 
    210 
     252                 objects += child.retrieve(get_vertices=get_vertices)
     253         return objects
     254       
    211255
    212256    def count(self, keywords=None):
  • anuga_core/source/anuga/utilities/test_quad.py

    r4654 r4655  
    11import unittest
     2from Numeric import array, allclose
     3
    24from quad import Cell, build_quadtree
    3 
    4 #from domain import *
    55from anuga.abstract_2d_finite_volumes.general_mesh import General_mesh as Mesh
    66
     
    2424
    2525        points = [a, b, c, d, e, f, g, h]
     26       
    2627        #bac, bce, ecf, dbe, daf, dae
    2728        vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]]
     
    5657        self.cell.split(4)
    5758
    58         result =  self.cell.search(x = 1, y = 101)
     59        result =  self.cell.search(x = 1, y = 101, get_vertices=True)
    5960        assert type(result) in [types.ListType,types.TupleType],\
    6061                                'should be a list'
     
    126127
    127128
    128         result = Q.search(3, 105)
     129        result = Q.search(3, 105, get_vertices=True)
    129130        assert type(result) in [types.ListType,types.TupleType],\
    130131                                'should be a list'
     
    151152        Q = build_quadtree(mesh)
    152153
     154    def test_retrieve_triangles(self):
     155
     156        cell = Cell(0, 6, 0, 6, 'cell', max_points_per_cell=4)
     157
     158        p0 = [2,1]
     159        p1 = [4,1]
     160        p2 = [4.,4]
     161        p3 = [2,4]
     162        p4 = [5,4]
     163
     164        points = [p0,p1,p2, p3, p4]
     165        #
     166        vertices = [[0,1,2],[0,2,3],[1,4,2]]
     167
     168        mesh = Mesh(points, vertices)
     169
     170        Q = build_quadtree(mesh)
     171        results = Q.search(5,1)
     172        assert len(results),2
     173        #print "results", results
     174        #print "results[0][0]", results[0][0]
     175        assert results[0],0
     176        assert results[1],2
     177        assert results[0][1],[[ 2.,  1.],
     178                     [ 4.,  1.],
     179                     [ 4.,  4.]]
     180        assert results[1][1],[[ 4.,  1.],
     181                     [ 5.,  4.],
     182                     [ 4.,  4.]]
     183        # this is the normals
     184        assert results[0][1][1],[[1.,  0.],
     185                     [-0.83205029,  0.5547002],
     186                     [ 0.,  -1.]]
     187                     
     188        # assert allclose(array(results),[[[ 2.,  1.],
     189        #[ 4.,  1.], [ 4.,  4.]], [[ 4.,  1.],[ 5.,  4.],[ 4.,  4.]]] )
     190        results = Q.search(5,4.)
     191        ### print "results",results
     192        # results_dic={}
     193        # results_dic.update(results)
     194        assert len(results),3
     195        #print "results_dic[0]", results_dic[0]
     196        assert results[0][1],[[ 2.,  1.],
     197                     [ 4.,  1.],
     198                     [ 4.,  4.]]
     199        assert results[1][1],[[ 2.,  1.],
     200                     [ 4.,  4.],
     201                     [ 2.,  4.]]
     202        assert results[2][1],[[ 4.,  1.],
     203                     [ 5.,  4.],
     204                     [ 4.,  4.]]
     205        #assert allclose(array(results),[[[ 2.,  1.],[ 4.,  1.], [ 4.,  4.]]
     206         #                               ,[[ 2.,  1.],[ 4.,  4.], [ 2.,  4.]],
     207        #[[ 4.,  1.],  [ 5.,  4.], [ 4.,  4.]],
     208         #                               [[ 4.,  1.], [ 5.,  4.], [ 4.,  4.]]])
     209       
     210       
    153211#-------------------------------------------------------------
    154212if __name__ == "__main__":
    155213
    156214    mysuite = unittest.makeSuite(Test_Quad,'test')
     215    # mysuite = unittest.makeSuite(Test_Quad,'test_retrieve_triangles')
    157216    runner = unittest.TextTestRunner()
    158217    runner.run(mysuite)
Note: See TracChangeset for help on using the changeset viewer.