Changeset 7718


Ignore:
Timestamp:
May 11, 2010, 8:12:13 PM (13 years ago)
Author:
hudson
Message:

Tag for new quadtree complete - beginning optimisations.

Location:
anuga_core/source/anuga
Files:
6 edited

Legend:

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

    r7317 r7718  
    2323import profile , pstats
    2424from math import sqrt
    25 
    26 from anuga.fit_interpolate.search_functions import search_times, \
    27      reset_search_times
    2825
    2926from anuga.fit_interpolate.interpolate import Interpolate
     
    262259        # backing up.
    263260       
    264         search_one_cell_time, search_more_cells_time = search_times()
    265         reset_search_times()
     261        #search_one_cell_time, search_more_cells_time = search_times()
     262        #reset_search_times()
    266263        #print "bench - search_one_cell_time",search_one_cell_time
    267264        #print "bench - search_more_cells_time", search_more_cells_time
  • anuga_core/source/anuga/fit_interpolate/test_search_functions.py

    r7716 r7718  
    33
    44import unittest
    5 from mesh_quadtree import MeshQuadtree, compute_interpolation_values
     5from mesh_quadtree import MeshQuadtree
    66
    77from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
     
    228228        assert k == 1
    229229       
    230        
    231     def test_compute_interpolation_values(self):
    232         """test_compute_interpolation_values
    233        
    234         Test that interpolation values are correc
    235        
    236         This test used to check element_found as output from
    237         find_triangle_compute_interpolation(triangle, n0, n1, n2, x)
    238         and that failed before 18th March 2009.
    239        
    240         Now this function no longer returns this flag, so the test
    241         is merely checknig the sigmas.
    242        
    243         """
    244        
    245         triangle = num.array([[306951.77151059, 6194462.14596986],
    246                               [306952.58403545, 6194459.65001246],
    247                               [306953.55109034, 6194462.0041216]])
    248 
    249 
    250         n0 = [0.92499377, -0.37998227]
    251         n1 = [0.07945684,  0.99683831]
    252         n2 = [-0.95088404, -0.30954732]
    253        
    254         x = [306953.344, 6194461.5]
    255        
    256         # Test that point is indeed inside triangle
    257         assert is_inside_polygon(x, triangle,
    258                                  closed=True, verbose=False)
    259         assert is_inside_triangle(x, triangle,
    260                                   closed=True, verbose=False)                                 
    261        
    262         sigma0, sigma1, sigma2 = \
    263             compute_interpolation_values(triangle, n0, n1, n2, x)
    264            
    265         msg = 'Point which is clearly inside triangle was not found'
    266         assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg
    267230
    268231#-------------------------------------------------------------
  • anuga_core/source/anuga/geometry/mesh_quadtree.py

    r7716 r7718  
    2727
    2828
    29 # FIXME(Ole): Could we come up with a less confusing structure?
    30 # FIXME(James): remove this global var
    3129LAST_TRIANGLE = [[-10, -10,
    3230                   (num.array([[max_float, max_float],
     
    5755        N = len(mesh)
    5856        self.mesh = mesh
    59         self.last_triangle = LAST_TRIANGLE            
     57        self.set_last_triangle() 
    6058
    6159        # Get x,y coordinates for all vertices for all triangles
     
    9189
    9290        """
     91       
     92        x = ensure_numeric(x, num.float)
     93                 
    9394        if self.last_triangle[0][1] != -10:
    9495            # check the last triangle found first
    9596            element_found, sigma0, sigma1, sigma2, k = \
    9697                       self._search_triangles_of_vertices(self.last_triangle, x)
    97 
    9898            if element_found:
    9999                return True, sigma0, sigma1, sigma2, k
    100100
    101101        branch = self.last_triangle[0][1]
    102        
    103         if branch == -10:
    104             branch = self   
    105102
    106103        # test neighbouring tris
     
    112109            return True, sigma0, sigma1, sigma2, k       
    113110
    114         # search to bottom of tree from last found leaf   
    115         tri_data = branch.search(x)
    116         triangles = self._trilist_from_data(tri_data)           
    117         element_found, sigma0, sigma1, sigma2, k = \
    118                     self._search_triangles_of_vertices(triangles, x)
    119         if element_found:
    120             return True, sigma0, sigma1, sigma2, k
    121 
    122111        # search rest of tree
    123112        element_found = False
    124         while branch:
    125             if not branch.parent:
    126                 # search from top of tree if we are at root
    127                 siblings = [self]
    128             else:
    129                 siblings = branch.get_siblings()
    130                
    131             for sibling in siblings:
     113        next_search = [branch]
     114        while branch:               
     115            for sibling in next_search:
    132116                tri_data = sibling.search(x)
    133117                triangles = self._trilist_from_data(tri_data)           
     
    136120                if element_found:
    137121                    return True, sigma0, sigma1, sigma2, k
    138                                        
     122           
     123            next_search = branch.get_siblings()                           
    139124            branch = branch.parent
    140125            if branch:
     
    144129                            self._search_triangles_of_vertices(triangles, x)
    145130                if element_found:
    146                     return True, sigma0, sigma1, sigma2, k       
     131                    return True, sigma0, sigma1, sigma2, k     
    147132
    148133        return element_found, sigma0, sigma1, sigma2, k
     
    158143        This function is responsible for most of the compute time in
    159144        fit and interpolate.
    160         """
    161         x = ensure_numeric(x, num.float)     
    162        
    163         # These statments are needed if triangles is empty
    164         sigma2 = -10.0
    165         sigma0 = -10.0
    166         sigma1 = -10.0
    167         k = -10
    168        
    169         # For all vertices in same cell as point x
    170         element_found = False   
     145        """ 
     146
    171147        for k, node, tri_verts_norms in triangles:
    172148            tri = tri_verts_norms[0]
     
    178154            # Input check disabled to speed things up.   
    179155            if bool(_is_inside_triangle(x, tri, int(True), 1.0e-12, 1.0e-12)):
    180                
    181156                n0, n1, n2 = tri_verts_norms[1]       
    182                 sigma0, sigma1, sigma2 =\
    183                     compute_interpolation_values(tri, n0, n1, n2, x)
    184                    
    185                 element_found = True   
     157                xi0, xi1, xi2 = tri
     158
     159                sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
     160                sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
     161                sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)             
    186162
    187163                # Don't look for any other triangles in the triangle list
    188164                self.last_triangle = [[k, node, tri_verts_norms]]
    189                 break           
    190                
    191         return element_found, sigma0, sigma1, sigma2, k
     165                return True, sigma0, sigma1, sigma2, k                     
     166        return False, -1, -1, -1, -10
    192167
    193168
    194169    def _trilist_from_data(self, indices):
    195         """return a list of lists. For the inner lists,
    196         The first element is the triangle index,
    197         the second element is a list.for this list
    198            the first element is a list of three (x, y) vertices,
    199            the following elements are the three triangle normals.
    200 
    201         """
    202 
    203170        ret_list = []
    204171        for i, node in indices:
     
    212179    def set_last_triangle(self):
    213180        self.last_triangle = LAST_TRIANGLE
     181        self.last_triangle[0][1] = self # point at root by default         
    214182
    215183   
    216 def compute_interpolation_values(triangle, n0, n1, n2, x):
    217     """Compute linear interpolation of point x and triangle.
    218    
    219     n0, n1, n2 are normal to the tree edges.
    220     """
    221184
    222     # Get the three vertex_points of candidate triangle k
    223     xi0, xi1, xi2 = triangle
    224 
    225     sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
    226     sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
    227     sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
    228 
    229     return sigma0, sigma1, sigma2
    230185               
  • anuga_core/source/anuga/geometry/quad.py

    r7714 r7718  
    1919    """
    2020 
    21     def __init__(self, extents, parent, depth = 0,
     21    def __init__(self, extents, parent,
    2222         name = 'cell'):
    2323 
     
    2727        self.extents = extents
    2828        self.parent = parent
    29         self.searched = [-10, -15]
    30         self.depth = depth
    3129       
    3230        # The points in this cell     
     
    3634   
    3735    def __repr__(self):
    38         str = '(%d)%s: leaves: %d' \
    39                % (self.depth, self.name , len(self.leaves))   
     36        str = '%s: leaves: %d' \
     37               % (self.name , len(self.leaves))   
    4038        if self.children:
    4139            str += ', children: %d' % (len(self.children))
     
    8482            # try splitting this cell and see if we get a trivial in
    8583            subregion1, subregion2 = self.extents.split()
    86             if subregion1.is_trivial_in(new_region):
    87                 self.children = [Cell(subregion1, self, depth=self.depth+1), Cell(subregion2, self, depth=self.depth+1)]   
    88                 self.children[0]._insert(new_leaf)
    89                 return
    90             elif subregion2.is_trivial_in(new_region):
    91                 self.children = [Cell(subregion1, self, depth=self.depth+1), Cell(subregion2, self, depth=self.depth+1)]   
    92                 self.children[1]._insert(new_leaf)
    93                 return               
     84           
     85            # option 1 - try splitting 4 ways
     86            subregion11, subregion12 = subregion1.split()                               
     87            subregion21, subregion22 = subregion2.split()
     88            regions = [subregion11, subregion12, subregion21, subregion22]
     89            for region in regions:
     90                if region.is_trivial_in(new_region):
     91                    self.children = [Cell(x, parent=self) for x in regions]
     92                    self._insert(new_leaf)
     93                    return               
     94
     95            # option 2 - try splitting 2 ways
     96            #if subregion1.is_trivial_in(new_region):
     97                #self.children = [Cell(subregion1, self), Cell(subregion2, self)]   
     98                #self.children[0]._insert(new_leaf)
     99                #return
     100            #elif subregion2.is_trivial_in(new_region):
     101                #self.children = [Cell(subregion1, self), Cell(subregion2, self)]   
     102                #self.children[1]._insert(new_leaf)
     103                #return               
    94104   
    95105        # recursion ended without finding a fit, so attach it as a leaf
     
    128138        if depth == 0:
    129139            log.critical()
    130         print '%s%s' % ('  '*depth, self.name), self.extents,' [', self.leaves, ']'
     140        print '%s%s'  % ('  '*depth, self.name), self.extents,' [', self.leaves, ']'
    131141        if self.children:
    132142            log.critical()
    133143            for child in self.children:
    134144                child.show(depth+1)
    135  
    136145
    137146    def search(self, x):
    138147        """return a list of possible intersections with geometry"""
    139        
     148     
    140149        intersecting_regions = self.test_leaves(x)
    141150       
  • anuga_core/source/anuga/geometry/test_meshquad.py

    r7716 r7718  
    145145        self.assertEqual(idx, 2)
    146146       
    147     def NOtest_get_parent():
    148         """ Get a child's parent.
     147    def NOtest_num_visits(self):
     148        """ Test optimisation code.
    149149        """
    150         cell = Cell(AABB(0, 6, 0, 6), 'cell')
     150        a = [3, 7]
     151        b = [5, 7]
     152        c = [5, 5]
     153        d = [7, 7]
     154        e = [15, 15]
     155        f = [15, 30]
     156        g = [30, 10]
     157        h = [30, 30]
    151158
    152         p0 = [2,1]
    153         p1 = [4,1]
    154         p2 = [4,4]
    155         p3 = [2,4]
    156         p4 = [5,4]
    157 
    158         points = [p0,p1,p2, p3, p4]
    159         vertices = [[0,1,2],[0,2,3],[1,4,2]]
     159        points = [a, b, c, d, e, f, g, h]
     160       
     161        #bac, bce, ecf, dbe, daf, dae
     162        vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]]
    160163
    161164        mesh = Mesh(points, vertices)
    162165
    163         Q = MeshQuadtree(mesh)
    164         results = Q.search([4.5, 3])
    165         assert len(results) == 1
    166         self.assertEqual(results[0], 2)
    167         results = Q.search([5,4.])
    168         self.assertEqual(len(results),1)
    169         self.assertEqual(results[0], 2)     
     166        Q = MeshQuadtree(mesh)   
     167
     168
     169        results = Q.search_fast([5.5, 5.5])
     170        print 'visits: ', Q.count_visits()
     171       
     172        Q.clear_visits()
     173        results = Q.search_fast([30, 10])
     174        print 'visits: ', Q.count_visits()
     175       
     176        print 'second time:'
     177
     178        Q.clear_visits()       
     179        results = Q.search_fast([5.5, 5.5])
     180        print 'visits: ', Q.count_visits()
    170181################################################################################
    171182
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r7712 r7718  
    18121812                gauge_values[i].append(w[i])
    18131813
    1814         #Reference (nautilus 26/6/2008)
    1815 
    1816         G0 = [-0.20000000000000001, -0.20000000000000001, -0.20000000000000001, -0.19729755264559332, -0.18143006781428656, -0.17264603126739803, -0.16490767093739186, -0.15580226478343825, -0.15341783464373857, -0.14986976789238646, -0.1483294354194323, -0.19385112169384708, -0.19907249891161938, -0.19911053885072424, -0.19913592111915288, -0.1991524040791994, -0.19916445999802629, -0.19917790655716, -0.19919280664859376, -0.19920527965237472, -0.19921381524768964, -0.1992178338453264, -0.19921637002152079, -0.19920783637167225, -0.19921088191840772, -0.19919037392917549, -0.19918812396862254, -0.19918939056535054, -0.19919336800395021, -0.19919832209289323, -0.19920329732763048, -0.19920824490004829, -0.19921301385072837, -0.19921748460975705, -0.1992215655528298, -0.19922516461276893, -0.19922820250665138, -0.19923065047412888, -0.19923254166903823, -0.1992339741831877, -0.19923510970454786, -0.19923612584283332, -0.19923716297385921, -0.199238301846201, -0.19923956605537632, -0.19924092368339749, -0.19924231511786381, -0.19924367609907154, -0.19924495328100339, -0.19924611576332413, -0.19924716014245455]
     1814        # Captured data from code manually inspected for correctness 11/5/2010
     1815
     1816        G0 = [-0.20000000000000001, -0.20000000000000001, -0.19999999639817381, -0.19900000000000001, -0.19787180620524364, -0.18775357859677891, -0.17717419377267535, -0.17038758570131016, -0.1632648895200452, -0.16084615293651969, -0.15796677906889614, -0.15610187077034218, -0.19890877268012952, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19767101056475031, -0.19664020763833073, -0.19625551309686509, -0.19622800460530576, -0.19670420796840654, -0.19757556964453959, -0.19853453054821305, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001, -0.19900000000000001]
    18171817        G1 = [-0.29999999999999993, -0.29999999999999993, -0.29312012778165086, -0.26736324399674027, -0.23794021958953263, -0.21878933349240881, -0.2063216250007715, -0.19630760483269205, -0.18778429435995911, -0.18241377240055098, -0.17860448861463396, -0.17352695896444634, -0.16449192186501546, -0.18953645096562932, -0.20557835913635514, -0.21071182844586539, -0.21191011610076083, -0.21068538990761979, -0.20737339293005266, -0.20407505794124087, -0.20190406545920248, -0.19987295968833538, -0.19832486432402194, -0.19774281454937118, -0.19763203737080062, -0.19769585527535841, -0.19795493497088903, -0.19841633252351859, -0.19886909990019369, -0.19940203310127189, -0.19990812587783294, -0.20016182994265733, -0.20026318580583596, -0.20022855664242351, -0.20010463117095048, -0.199930824934102, -0.19973147245180795, -0.19957882377401115, -0.19948188719957555, -0.19944190810764481, -0.19945196574982568, -0.19949690579414533, -0.19955788251561279, -0.19961580710612042, -0.19965570581001513, -0.19966924072268322, -0.19966164565962002, -0.1996429044375452, -0.19962419292036049, -0.19961203640143024, -0.19960927744740292]
    18181818        G2 = [-0.40000000000000002, -0.39653752812451043, -0.32890250043784469, -0.29967800724628818, -0.27572700190320054, -0.25893877811152488, -0.2422572785756886, -0.22758771043596362, -0.21433402198450011, -0.20266876405969875, -0.19393588052745447, -0.18766846928826031, -0.18182722392707523, -0.16865402349317951, -0.18110962350013771, -0.19547215504627991, -0.20229670797148233, -0.20494160396870145, -0.20614855426298931, -0.20616358034144364, -0.2050749944506553, -0.20329491765401383, -0.20174840479413605, -0.20040155802388543, -0.19945041001985703, -0.19898209259951644, -0.19868944730468457, -0.19849574784901269, -0.19855116942355155, -0.19885365879108613, -0.19925381355621719, -0.19971339601933866, -0.20000203652933318, -0.20011337898513715, -0.2001188815225943, -0.20007116187977508, -0.20000646342088527, -0.19993034594397555, -0.19985512250322437, -0.19978530424723981, -0.19973800656091345, -0.19971854873923284, -0.19972327598147538, -0.19974189986453444, -0.19976101408429672, -0.19977705650936498, -0.19978787136068635, -0.19979206946039513, -0.19978865745584601, -0.19978181317853855, -0.19977875136946421]
Note: See TracChangeset for help on using the changeset viewer.