Changeset 608


Ignore:
Timestamp:
Nov 22, 2004, 3:59:43 PM (20 years ago)
Author:
duncan
Message:

change least_squares algorithm, so it increases its cell search, if the point is not within any of the initial cells triangles. 1st cut.

Location:
inundation/ga/storm_surge/pyvolution
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/interpolate_sww.py

    r600 r608  
    7373        self.point_coordinates = self.read_xya(file_name)
    7474        self.interp.build_interpolation_matrix_A(self.point_coordinates)
    75         #print "self.interp.A",self.interp.A
    7675
    7776        #Change this if you want to generalise this application -
  • inundation/ga/storm_surge/pyvolution/least_squares.py

    r606 r608  
    239239            #Find vertices near x
    240240            candidate_vertices = root.search(x[0], x[1])
    241 
     241               
     242            element_found, sigma0, sigma1, sigma2, k = \
     243                self.search_triangles_of_vertices(candidate_vertices, x)
     244
     245            is_more_elements = True
     246            while not element_found and is_more_elements:
     247                candidate_vertices = root.expand_search()
     248                if candidate_vertices == []:
     249                    # All the triangles have been searched.
     250                    is_more_elements = False
     251                else:
     252                    element_found, sigma0, sigma1, sigma2, k = \
     253                      self.search_triangles_of_vertices(candidate_vertices, x)
     254                   
     255               
    242256           
     257            #Update interpolation matrix A if necessary     
     258            if element_found is True:       
     259                #Assign values to matrix A
     260
     261                j0 = self.mesh.triangles[k,0] #Global vertex id
     262                #self.A[i, j0] = sigma0
     263
     264                j1 = self.mesh.triangles[k,1] #Global vertex id
     265                #self.A[i, j1] = sigma1
     266
     267                j2 = self.mesh.triangles[k,2] #Global vertex id
     268                #self.A[i, j2] = sigma2
     269
     270                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
     271                js     = [j0,j1,j2]
     272
     273                for j in js:
     274                    self.A[i,j] = sigmas[j]
     275                    for k in js:
     276                        self.AtA[j,k] += sigmas[j]*sigmas[k]
     277            else:
     278                pass
     279                #Ok if there is no triangle for datapoint
     280                #(as in brute force version)
     281                #raise 'Could not find triangle for point', x
     282
     283
     284    def search_triangles_of_vertices(self, candidate_vertices, x):
    243285            #Find triangle containing x:
    244             element_found = False           
    245            
     286            element_found = False
     287
     288            # This will be returned if element_found = False
     289            sigma2 = -10.0
     290            sigma0 = -10.0
     291            sigma1 = -10.0
     292
    246293            #For all vertices in same cell as point x
    247294            for v in candidate_vertices:
     
    255302                    xi2 = self.mesh.get_vertex_coordinate(k, 2)     
    256303
     304                    #print "PDSG - k", k
     305                    #print "PDSG - xi0", xi0
     306                    #print "PDSG - xi1", xi1
     307                    #print "PDSG - xi2", xi2
     308                    #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))" \
     309                    #      % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
     310                   
    257311                    #Get the three normals
    258312                    n0 = self.mesh.get_normal(k, 0)
     
    267321                    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
    268322
     323                    #print "PDSG - sigma0", sigma0
     324                    #print "PDSG - sigma1", sigma1
     325                    #print "PDSG - sigma2", sigma2
     326                   
    269327                    #FIXME: Maybe move out to test or something
    270328                    epsilon = 1.0e-6
     
    275333                    #Sigmas can get negative within
    276334                    #machine precision on some machines (e.g nautilus)
    277                     #Hence the small eps
    278                    
     335                    #Hence the small eps                   
    279336                    eps = 1.0e-15
    280337                    if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
     
    285342                    #Don't look for any other triangle
    286343                    break
     344            return element_found, sigma0, sigma1, sigma2, k     
    287345                   
    288 
    289 
    290            
    291             #Update interpolation matrix A if necessary     
    292             if element_found is True:       
    293                 #Assign values to matrix A
    294 
    295                 j0 = self.mesh.triangles[k,0] #Global vertex id
    296                 #self.A[i, j0] = sigma0
    297 
    298                 j1 = self.mesh.triangles[k,1] #Global vertex id
    299                 #self.A[i, j1] = sigma1
    300 
    301                 j2 = self.mesh.triangles[k,2] #Global vertex id
    302                 #self.A[i, j2] = sigma2
    303 
    304                 sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
    305                 js     = [j0,j1,j2]
    306 
    307                 for j in js:
    308                     self.A[i,j] = sigmas[j]
    309                     for k in js:
    310                         self.AtA[j,k] += sigmas[j]*sigmas[k]
    311             else:
    312                 pass
    313                 #Ok if there is no triangle for datapoint
    314                 #(as in brute force version)
    315                 #raise 'Could not find triangle for point', x
    316 
    317346
    318347       
  • inundation/ga/storm_surge/pyvolution/quad.py

    r484 r608  
    7070        self.AddChild(Cell(cs,(cn+cs)/2,(cw+ce)/2,ce,self.name+'_se'))
    7171        self.AddChild(Cell(cs,(cn+cs)/2,cw,(cw+ce)/2,self.name+'_sw'))
    72      
     72        
    7373 
    74 
    7574    def search(self, x, y):
     75    #def search_new(self, x, y):
    7676        """Find all point indices sharing the same cell as point (x, y)
    7777        """
    78      
     78        branch = []
    7979        points = []
    8080        if self.children:
    8181            for child in self:
    8282                if child.contains(x,y):
    83                      points += child.search(x,y)
     83                    branch.append(self)
     84                    points, branch = child.search_branch(x,y, branch)
    8485        else:
    8586            # Leaf node: Get actual waypoints
    8687            points = self.retrieve()
    87            
     88
     89        self.branch = branch   
    8890        return points
    8991
     92
     93    def search_branch(self, x, y, branch):
     94        """Find all point indices sharing the same cell as point (x, y)
     95        """
     96        points = []
     97        if self.children:
     98            for child in self:
     99                if child.contains(x,y):
     100                    branch.append(self)
     101                    points, branch = child.search_branch(x,y, branch)
     102                   
     103        else:
     104            # Leaf node: Get actual waypoints
     105            points = self.retrieve()     
     106        return points, branch
     107
     108
     109    def expand_search(self):
     110        """Find all point indices 'up' one cell from the last search
     111        """
     112        points = []
     113        if self.branch == []:
     114            points = []
     115        else:
     116            larger_cell = self.branch.pop()
     117            #print "larger_cell ", larger_cell.show() # Get_tree()
     118            for child in larger_cell:
     119                #This means it will also go down the branch known to be bad.
     120                #So, if things are too slow, this could be speed up
     121                points += child.retrieve()
     122        return points
    90123
    91124
     
    367400    #Hence, the root cell needs to be expanded slightly
    368401    ymax += (ymax-ymin)/10
    369     xmax += (xmax-xmin)/10    
     402    xmax += (xmax-xmin)/10   
    370403
    371404
  • inundation/ga/storm_surge/pyvolution/test_least_squares.py

    r570 r608  
    4040       
    4141
     42    def test_quad_tree(self):
     43        p0 = [-10.0, -10.0]
     44        p1 = [20.0, -10.0]
     45        p2 = [-10.0, 20.0]
     46        p3 = [10.0, 50.0]
     47        p4 = [30.0, 30.0]
     48        p5 = [50.0, 10.0]
     49        p6 = [40.0, 60.0]
     50        p7 = [60.0, 40.0]
     51        p8 = [-66.0, 20.0]
     52        p9 = [10.0, -66.0]
     53       
     54        points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
     55        triangles = [ [0, 1, 2],
     56                      [3, 2, 4],
     57                      [4, 2, 1],
     58                      [4, 1, 5],
     59                      [3, 4, 6],
     60                      [6, 4, 7],
     61                      [7, 4, 5],
     62                      [8, 0, 2],
     63                      [0, 9, 1]]   
     64
     65        data = [ [4,4] ]
     66       
     67        interp = Interpolation(points, triangles, data, alpha = 0.0)
     68        #print "PDSG - interp.get_A()", interp.get_A()
     69        answer =  [ [ 0.06666667,  0.46666667,  0.46666667,  0.,
     70                      0., 0. , 0., 0., 0., 0.]]
     71        assert allclose(interp.get_A(), answer)
     72       
    4273
    4374    def test_datapoints_at_vertices(self):
Note: See TracChangeset for help on using the changeset viewer.