Changeset 2190

Timestamp:

Jan 9, 2006, 5:17:32 PM (19 years ago)

Author:

duncan

Message:

refactoring out search functions

Location:

inundation/fit_interpolate

Files:

• interpolate.py (modified) (5 diffs)
• search_functions.py (added)
• test_interpolate.py (modified) (1 diff)

inundation/fit_interpolate/interpolate.py

@@ -22 +22 @@
 from pyvolution.cg_solve import conjugate_gradient, VectorShapeError
 from coordinate_transforms.geo_reference import Geo_reference
-
 from pyvolution.quad import build_quadtree
 from utilities.numerical_tools import ensure_numeric
 from utilities.polygon import inside_polygon
 
+from search_functions import search_tree_of_vertices
 
 class Interpolate:
@@ -138 +138 @@
         A = Sparse(n,m)
 
-
-        # I think this (along with other expanded_quad_searches stuff) can be removed 
-        #self.expanded_quad_searches = []
-
         #Compute matrix elements
         for i in range(n):
@@ -151 +147 @@
             
             element_found, sigma0, sigma1, sigma2, k = \
-                           self._search_tree_of_vertices(x)    
+                           search_tree_of_vertices(self.root, self.mesh, x)
+            
             #Update interpolation matrix A if necessary
             if element_found is True:
@@ -170 +167 @@
         return A
 
-    def _search_tree_of_vertices(self, x):
+    def _search_tree_of_vertices(self, root, mesh, x):
         """
         Find the triangle (element) that the point x is in.
 
+        root: A quad tree of the vertices
         Return the associated sigma and k values
            (and if the element was found) .
@@ -187 +185 @@
             
         #Find vertices near x
-        candidate_vertices = self.root.search(x[0], x[1])
+        candidate_vertices = root.search(x[0], x[1])
         is_more_elements = True
 
         element_found, sigma0, sigma1, sigma2, k = \
-                       self._search_triangles_of_vertices(candidate_vertices, x)
-        #FIXME Do we need this var?
-        # Exclude points outside the mesh before removing this. 
-        #first_expansion = True
+                       self._search_triangles_of_vertices(mesh,
+                                                          candidate_vertices, x)
         while not element_found and is_more_elements:
-            #if verbose: print 'Expanding search'
-            #if first_expansion == True:
-            #self.expanded_quad_searches.append(1)
-            #first_expansion = False
-            #else:
-            #end = len(self.expanded_quad_searches) - 1
-            #assert end >= 0
-            #self.expanded_quad_searches[end] += 1
-            candidate_vertices, branch = self.root.expand_search()
+            candidate_vertices, branch = root.expand_search()
             if branch == []:
                 # Searching all the verts from the root cell that haven't
                 # been searched.  This is the last try
                 element_found, sigma0, sigma1, sigma2, k = \
-                      self._search_triangles_of_vertices(candidate_vertices, x)
+                      self._search_triangles_of_vertices(mesh,
+                                                         candidate_vertices, x)
                 is_more_elements = False
             else:
                 element_found, sigma0, sigma1, sigma2, k = \
-                      self._search_triangles_of_vertices(candidate_vertices, x)
+                      self._search_triangles_of_vertices(mesh,
+                                                         candidate_vertices, x)
 
         return element_found, sigma0, sigma1, sigma2, k
 
-    def _search_triangles_of_vertices(self,
-                                 candidate_vertices, x):
-            #Find triangle containing x:
-            element_found = False
-
-            # This will be returned if element_found = False
-            sigma2 = -10.0
-            sigma0 = -10.0
-            sigma1 = -10.0
-            k = -10.0
-            #print "*$* candidate_vertices", candidate_vertices
-            #For all vertices in same cell as point x
-            for v in candidate_vertices:
-                #FIXME (DSG-DSG): this catches verts with no triangle.
-                #Currently pmesh is producing these.
-                #this should be stopped, 
-                if self.mesh.vertexlist[v] is None:
-                    continue
-                #for each triangle id (k) which has v as a vertex
-                for k, _ in self.mesh.vertexlist[v]:
-                    #Get the three vertex_points of candidate triangle
-                    xi0 = self.mesh.get_vertex_coordinate(k, 0)
-                    xi1 = self.mesh.get_vertex_coordinate(k, 1)
-                    xi2 = self.mesh.get_vertex_coordinate(k, 2)
-
-                    #Get the three normals
-                    n0 = self.mesh.get_normal(k, 0)
-                    n1 = self.mesh.get_normal(k, 1)
-                    n2 = self.mesh.get_normal(k, 2)
-
-                    #Compute interpolation
-                    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
-                    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
-                    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
-
-                    #FIXME: Maybe move out to test or something
-                    epsilon = 1.0e-6
-                    assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
-
-                    #Check that this triangle contains the data point
-
-                    #Sigmas can get negative within
-                    #machine precision on some machines (e.g nautilus)
-                    #Hence the small eps
-                    eps = 1.0e-15
-                    if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
-                        element_found = True
-                        break
-
-                if element_found is True:
-                    #Don't look for any other triangle
+    def _search_triangles_of_vertices(self, mesh, candidate_vertices, x):
+        #Find triangle containing x:
+        element_found = False
+
+        # This will be returned if element_found = False
+        sigma2 = -10.0
+        sigma0 = -10.0
+        sigma1 = -10.0
+        k = -10.0
+        #print "*$* candidate_vertices", candidate_vertices
+        #For all vertices in same cell as point x
+        for v in candidate_vertices:
+            #FIXME (DSG-DSG): this catches verts with no triangle.
+            #Currently pmesh is producing these.
+            #this should be stopped, 
+            if mesh.vertexlist[v] is None:
+                continue
+            #for each triangle id (k) which has v as a vertex
+            for k, _ in mesh.vertexlist[v]:
+                #Get the three vertex_points of candidate triangle
+                xi0 = mesh.get_vertex_coordinate(k, 0)
+                xi1 = mesh.get_vertex_coordinate(k, 1)
+                xi2 = mesh.get_vertex_coordinate(k, 2)
+
+                #Get the three normals
+                n0 = mesh.get_normal(k, 0)
+                n1 = mesh.get_normal(k, 1)
+                n2 = mesh.get_normal(k, 2)
+
+                #Compute interpolation
+                sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
+                sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
+                sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
+                    
+                #FIXME: Maybe move out to test or something
+                epsilon = 1.0e-6
+                assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
+
+                #Check that this triangle contains the data point
+                
+                #Sigmas can get negative within
+                #machine precision on some machines (e.g nautilus)
+                #Hence the small eps
+                eps = 1.0e-15
+                if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
+                    element_found = True
                     break
-            return element_found, sigma0, sigma1, sigma2, k
+
+            if element_found is True:
+                #Don't look for any other triangle
+                break
+        return element_found, sigma0, sigma1, sigma2, k
 
 

inundation/fit_interpolate/test_interpolate.py

@@ -235 +235 @@
 
 
-        results = interp._build_interpolation_matrix_A(data).todense()
-        assert allclose(results, answer)
+        #results = interp._build_interpolation_matrix_A(data).todense()
+
+        assert allclose(A, answer)