Changeset 4654


Ignore:
Timestamp:
Aug 6, 2007, 10:02:39 AM (18 years ago)
Author:
duncan
Message:

did "svn merge -r 4653:4652 ​https://datamining.anu.edu.au/svn/ga" in the head dir, to remove the optimisation until it has been benchmarked more.

Location:
anuga_core/source/anuga
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r4653 r4654  
    3737    k = -10.0
    3838           
    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])
     39    #Find vertices near x
     40    candidate_vertices = root.search(x[0], x[1])
    4341    is_more_elements = True
    4442
    4543    element_found, sigma0, sigma1, sigma2, k = \
    4644                   _search_triangles_of_vertices(mesh,
    47                                                  triangles, x)
     45                                                 candidate_vertices, x)
    4846    while not element_found and is_more_elements:
    49         triangles, branch = root.expand_search()
     47        candidate_vertices, branch = root.expand_search()
    5048        if branch == []:
    5149            # Searching all the verts from the root cell that haven't
    5250            # been searched.  This is the last try
    5351            element_found, sigma0, sigma1, sigma2, k = \
    54                            _search_triangles_of_vertices(mesh,triangles, x)
     52                           _search_triangles_of_vertices(mesh,
     53                                                         candidate_vertices, x)
    5554            is_more_elements = False
    5655        else:
    5756            element_found, sigma0, sigma1, sigma2, k = \
    58                        _search_triangles_of_vertices(mesh,triangles, x)
     57                       _search_triangles_of_vertices(mesh,
     58                                                     candidate_vertices, x)
    5959
    6060    return element_found, sigma0, sigma1, sigma2, k
    6161
    6262
    63 def _search_triangles_of_vertices(mesh, triangles, x):
     63def _search_triangles_of_vertices(mesh, candidate_vertices, 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 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
     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
    8998        # Find triangle that contains x (if any) and interpolate
    90         element_found, sigma0, sigma1, sigma2 =\
    91                        find_triangle_compute_interpolation(tri, n0, n1, n2, x)
     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
    92110        if element_found is True:
    93             # Don't look for any other triangles in the triangle list
     111            # Don't look for any other triangle_lists from the
     112            # candidate_vertices
    94113            break
     114       
    95115    return element_found, sigma0, sigma1, sigma2, k
    96116
    97117
    98118           
    99 def find_triangle_compute_interpolation(triangle, n0, n1, n2, x):
     119def find_triangle_compute_interpolation(mesh, k, x):
    100120    """Compute linear interpolation of point x and triangle k in mesh.
    101121    It is assumed that x belongs to triangle k.
    … …  
    103123
    104124    # Get the three vertex_points of candidate triangle k
    105     xi0, xi1, xi2 = triangle
     125    xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)
    106126
    107127    # this is where we can call some fast c code.
    … …  
    118138   
    119139    if  x[0] > xmax + epsilon:
    120         return False,0,0,0
     140        return False,0,0,0,0
    121141    if  x[0] < xmin - epsilon:
    122         return False,0,0,0
     142        return False,0,0,0,0
    123143    if  x[1] > ymax + epsilon:
    124         return False,0,0,0
     144        return False,0,0,0,0
    125145    if  x[1] < ymin - epsilon:
    126         return False,0,0,0
     146        return False,0,0,0,0
    127147   
    128148    # Get the three normals
    129     #n0 = norms[0] 
    130     #n1 = norms[1]
    131     #n2 = norms[2]
     149    n0 = mesh.get_normal(k, 0)
     150    n1 = mesh.get_normal(k, 1)
     151    n2 = mesh.get_normal(k, 2)           
     152       
    132153       
    133154    # Compute interpolation
    … …  
    150171    else:
    151172        element_found = False
    152     return element_found, sigma0, sigma1, sigma2
     173    return element_found, sigma0, sigma1, sigma2, k

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

    r4653 r4654  
    1414ofile = 'lbm_resultsII.csv'
    1515delimiter = ','
    16 run_profile = False #True
     16run_profile = True #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

    r4653 r4654  
    7272       
    7373 
    74     def search(self, x, y, get_vertices=False):
     74    def search(self, x, y):
     75    #def search_new(self, x, y):
    7576        """Find all point indices sharing the same cell as point (x, y)
    7677        """
    … …  
    8384                    brothers.remove(child)
    8485                    branch.append(brothers)
    85                     points, branch = child.search_branch(x,y, branch,
    86                                                   get_vertices=get_vertices)
     86                    points, branch = child.search_branch(x,y, branch)
    8787        else:
    8888            # Leaf node: Get actual waypoints
    89             points = self.retrieve(get_vertices=get_vertices)
     89            points = self.retrieve()
     90
    9091        self.branch = branch   
    9192        return points
    9293
    9394
    94     def search_branch(self, x, y, branch, get_vertices=False):
     95    def search_branch(self, x, y, branch):
    9596        """Find all point indices sharing the same cell as point (x, y)
    9697        """
    … …  
    102103                    brothers.remove(child)
    103104                    branch.append(brothers)
    104                     points, branch = child.search_branch(x,y, branch,
    105                                                   get_vertices=get_vertices)
     105                    points, branch = child.search_branch(x,y, branch)
    106106                   
    107107        else:
    108108            # Leaf node: Get actual waypoints
    109             points = self.retrieve(get_vertices=get_vertices)     
     109            points = self.retrieve()     
    110110        return points, branch
    111111
    112112
    113     def expand_search(self, get_vertices=False):
     113    def expand_search(self):
    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(get_vertices=get_vertices)
     123                points += cell.retrieve()
    124124        return points, self.branch
    125125
    … …  
    200200
    201201
    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):
     202    def retrieve(self):
    233203         objects = []
    234204         if self.children is None:
    … …  
    239209         return objects 
    240210
    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        
    248211
    249212    def count(self, keywords=None):

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

    r4653 r4654  
    11import unittest
    2 from Numeric import array, allclose
     2from quad import Cell, build_quadtree
    33
    4 from quad import Cell, build_quadtree
     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        
    2726        #bac, bce, ecf, dbe, daf, dae
    2827        vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]]
    … …  
    5756        self.cell.split(4)
    5857
    59         result =  self.cell.search(x = 1, y = 101, get_vertices=True)
     58        result =  self.cell.search(x = 1, y = 101)
    6059        assert type(result) in [types.ListType,types.TupleType],\
    6160                                'should be a list'
    … …  
    127126
    128127
    129         result = Q.search(3, 105, get_vertices=True)
     128        result = Q.search(3, 105)
    130129        assert type(result) in [types.ListType,types.TupleType],\
    131130                                'should be a list'
    … …  
    152151        Q = build_quadtree(mesh)
    153152
    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        
    210153#-------------------------------------------------------------
    211154if __name__ == "__main__":
    212155
    213156    mysuite = unittest.makeSuite(Test_Quad,'test')
    214     # mysuite = unittest.makeSuite(Test_Quad,'test_retrieve_triangles')
    215157    runner = unittest.TextTestRunner()
    216158    runner.run(mysuite)
Note: See TracChangeset for help on using the changeset viewer.