Changeset 4653


Ignore:
Timestamp:
Aug 6, 2007, 9:56:50 AM (17 years ago)
Author:
duncan
Message:

checking in for benchmarking. 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 thing up.

Location:
anuga_core/source/anuga
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r4651 r4653  
    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/ticket178_benchmark.py

    r4652 r4653  
    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

    r3945 r4653  
    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        """
    … …  
    121121            for cell in three_cells:
    122122                #print "cell ", cell.show()
    123                 points += cell.retrieve()
     123                points += cell.retrieve(get_vertices=get_vertices)
    124124        return points, self.branch
    125125
    … …  
    200200
    201201
    202     def retrieve(self):
     202    def retrieve_triangles(self):
     203        """return a list of lists. For the inner lists,
     204        The first element is the triangle index,
     205        the second element is a list.for this list
     206           the first element is a list of three (x, y) vertices,
     207           the following elements are the three triangle normals.
     208
     209        This info is used in searching for a triangle that a point is in.
     210        """
     211        # FIXME Tidy up the structure that is returned.
     212        # if the triangles att has been made
     213        # return it.
     214        if not hasattr(self,'triangles'):
     215            # use a dictionary to remove duplicates
     216            triangles = {}
     217            verts = self.retrieve_vertices()
     218            # print "verts", verts
     219            for vert in verts:
     220                triangle_list = self.__class__.mesh.get_triangles_and_vertices_per_node(vert)
     221                for k, _ in triangle_list:
     222                    if not triangles.has_key(k):
     223                        # print 'k',k
     224                        tri = self.__class__.mesh.get_vertex_coordinates(k)
     225                        n0 = self.__class__.mesh.get_normal(k, 0)
     226                        n1 = self.__class__.mesh.get_normal(k, 1)
     227                        n2 = self.__class__.mesh.get_normal(k, 2)
     228                        triangles[k]=(tri, (n0, n1, n2))
     229            self.triangles = triangles.items()
     230        return self.triangles
     231           
     232    def retrieve_vertices(self):
    203233         objects = []
    204234         if self.children is None:
    … …  
    209239         return objects 
    210240
     241
     242    def retrieve(self, get_vertices=True):
     243        if get_vertices is True:
     244            return self.retrieve_vertices()
     245        else:
     246            return self.retrieve_triangles()
     247       
    211248
    212249    def count(self, keywords=None):

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

    r3945 r4653  
    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       
    153210#-------------------------------------------------------------
    154211if __name__ == "__main__":
    155212
    156213    mysuite = unittest.makeSuite(Test_Quad,'test')
     214    # mysuite = unittest.makeSuite(Test_Quad,'test_retrieve_triangles')
    157215    runner = unittest.TextTestRunner()
    158216    runner.run(mysuite)
Note: See TracChangeset for help on using the changeset viewer.