Changeset 4655

Timestamp:

Aug 6, 2007, 1:53:26 PM (18 years ago)

Author:

duncan

Message:

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 things up. Also, the old point info is deleted to save memory.

Location:

anuga_core/source/anuga

Files:

• fit_interpolate/search_functions.py (modified) (5 diffs)
• fit_interpolate/test_search_functions.py (modified) (2 diffs)
• fit_interpolate/ticket178_benchmark.py (modified) (1 diff)
• utilities/quad.py (modified) (5 diffs)
• utilities/test_quad.py (modified) (5 diffs)

anuga_core/source/anuga/fit_interpolate/search_functions.py

@@ -37 +37 @@
     k = -10.0
             
-    #Find vertices near x
-    candidate_vertices = root.search(x[0], x[1])
+    #Get triangles in the cell that the point is in.
+    # Triangle is a list, first element triangle_id,
+    # second element the triangle
+    triangles = root.search(x[0], x[1])
     is_more_elements = True
 
     element_found, sigma0, sigma1, sigma2, k = \
                    _search_triangles_of_vertices(mesh,
-                                                 candidate_vertices, x)
+                                                 triangles, x)
     while not element_found and is_more_elements:
-        candidate_vertices, branch = root.expand_search()
+        triangles, 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 = \
-                           _search_triangles_of_vertices(mesh,
-                                                         candidate_vertices, x)
+                           _search_triangles_of_vertices(mesh,triangles, x)
             is_more_elements = False
         else:
             element_found, sigma0, sigma1, sigma2, k = \
-                       _search_triangles_of_vertices(mesh,
-                                                     candidate_vertices, x)
+                       _search_triangles_of_vertices(mesh,triangles, x)
 
     return element_found, sigma0, sigma1, sigma2, k
 
 
-def _search_triangles_of_vertices(mesh, candidate_vertices, x):
+def _search_triangles_of_vertices(mesh, triangles, x):
     """Search for triangle containing x amongs candidate_vertices in mesh
 
@@ -82 +82 @@
     
     #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.number_of_triangles_per_node[v] == 0:
-            continue
-        
-        # Get all triangles which has v as a vertex
-        # The list has elements (triangle, vertex), but only the
-        # first component will be used here
-        triangle_list = mesh.get_triangles_and_vertices_per_node(node=v)
-
+    for k, tri_verts_norms in triangles:
+        tri = tri_verts_norms[0]
+        n0, n1, n2 = tri_verts_norms[1]
+        # k is the triangle index
+        # tri is a list of verts (x, y), representing a tringle
         # Find triangle that contains x (if any) and interpolate
-        
-        for k, _ in triangle_list:
-            element_found, sigma0, sigma1, sigma2, k =\
-                           find_triangle_compute_interpolation(mesh,
-                                                               k,
-                                                               x)
-
-            if element_found is True:
-                # Don't look for any other triangles in the triangle list
-                break
-
+        element_found, sigma0, sigma1, sigma2 =\
+                       find_triangle_compute_interpolation(tri, n0, n1, n2, x)
         if element_found is True:
-            # Don't look for any other triangle_lists from the
-            # candidate_vertices
+            # Don't look for any other triangles in the triangle list
             break
-        
     return element_found, sigma0, sigma1, sigma2, k
 
 
             
-def find_triangle_compute_interpolation(mesh, k, x):
+def find_triangle_compute_interpolation(triangle, n0, n1, n2, x):
     """Compute linear interpolation of point x and triangle k in mesh.
     It is assumed that x belongs to triangle k.
@@ -123 +103 @@
 
     # Get the three vertex_points of candidate triangle k
-    xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)
+    xi0, xi1, xi2 = triangle
 
     # this is where we can call some fast c code.
@@ -138 +118 @@
     
     if  x[0] > xmax + epsilon:
-        return False,0,0,0,0
+        return False,0,0,0
     if  x[0] < xmin - epsilon:
-        return False,0,0,0,0
+        return False,0,0,0
     if  x[1] > ymax + epsilon:
-        return False,0,0,0,0
+        return False,0,0,0
     if  x[1] < ymin - epsilon:
-        return False,0,0,0,0
+        return False,0,0,0
     
     # Get the three normals 
-    n0 = mesh.get_normal(k, 0)
-    n1 = mesh.get_normal(k, 1)
-    n2 = mesh.get_normal(k, 2)            
-        
+    #n0 = norms[0]  
+    #n1 = norms[1]
+    #n2 = norms[2]
         
     # Compute interpolation
@@ -171 +150 @@
     else:
         element_found = False 
-    return element_found, sigma0, sigma1, sigma2, k
+    return element_found, sigma0, sigma1, sigma2

anuga_core/source/anuga/fit_interpolate/test_search_functions.py

@@ -12 +12 @@
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
 from anuga.utilities.polygon import is_inside_polygon
-from anuga.utilities.quad import build_quadtree
-
-
-
+from anuga.utilities.quad import build_quadtree, Cell
 
 
@@ -161 +158 @@
                 assert found is True
             else:
-                assert found is False   
+                assert found is False
+
+                
+
+    def expanding_search(self):
+        """test_larger mesh and different quad trees
+        """
+        
+        p0 = [2,1]
+        p1 = [4,1]
+        p2 = [4.,4]
+        p3 = [2,4]
+        p4 = [5,4]
+
+        p5 = [-1,-1]
+        p6 = [1,-1]
+        p7 = [1,1]
+        p8 = [-1,1]
+
+        points = [p0,p1,p2, p3,p4,p5,p6,p7,p8]
+        #
+        vertices = [[0,1,2],[0,2,3],[1,4,2],[5,6,7], [5,7,8]]
+        mesh = Mesh(points, vertices)
+
+        # Don't do this, want to control the mx and mins
+        #root = build_quadtree(mesh, max_points_per_cell=4)
+    
+        #Initialise
+        Cell.initialise(mesh)
+
+        root = Cell(-3, 9, -3, 9,
+                    max_points_per_cell = 4)
+        #Insert indices of all vertices
+        root.insert( range(mesh.number_of_nodes) )
+
+        #Build quad tree and return
+        root.split()
+        
+        # One point
+        #x = [3.5, 1.5]
+        x = [2.5, 1.5]
+        element_found, sigma0, sigma1, sigma2, k = \
+                       search_tree_of_vertices(root, mesh, x)
+        # One point
+        x = [3.00005, 2.999994]
+        element_found, sigma0, sigma1, sigma2, k = \
+                       search_tree_of_vertices(root, mesh, x)
+        assert element_found is True
+        assert k == 1
+        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_search_functions,'test')
+    #suite = unittest.makeSuite(Test_search_functions,'expanding_search')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)

anuga_core/source/anuga/fit_interpolate/ticket178_benchmark.py

@@ -14 +14 @@
 ofile = 'lbm_resultsII.csv'
 delimiter = ','
-run_profile = True #False #True
+run_profile = False #True
 is_fit_list = [True, False]
 num_of_points_list = [3, 200, 600, 2000, 6000, 10000, 20000] 

anuga_core/source/anuga/utilities/quad.py

@@ -72 +72 @@
         
  
-    def search(self, x, y):
-    #def search_new(self, x, y):
+    def search(self, x, y, get_vertices=False):
         """Find all point indices sharing the same cell as point (x, y)
         """
@@ -84 +83 @@
                     brothers.remove(child)
                     branch.append(brothers)
-                    points, branch = child.search_branch(x,y, branch)
+                    points, branch = child.search_branch(x,y, branch,
+                                                  get_vertices=get_vertices)
         else:
             # Leaf node: Get actual waypoints
-            points = self.retrieve()
-
+            points = self.retrieve(get_vertices=get_vertices)
         self.branch = branch   
         return points
 
 
-    def search_branch(self, x, y, branch):
+    def search_branch(self, x, y, branch, get_vertices=False):
         """Find all point indices sharing the same cell as point (x, y)
         """
@@ -103 +102 @@
                     brothers.remove(child)
                     branch.append(brothers)
-                    points, branch = child.search_branch(x,y, branch)
+                    points, branch = child.search_branch(x,y, branch,
+                                                  get_vertices=get_vertices)
                     
         else:
             # Leaf node: Get actual waypoints
-            points = self.retrieve()      
+            points = self.retrieve(get_vertices=get_vertices)      
         return points, branch
 
 
-    def expand_search(self):
+    def expand_search(self, get_vertices=False):
         """Find all point indices 'up' one cell from the last search
         """
+        
         points = []
         if self.branch == []:
@@ -121 +122 @@
             for cell in three_cells:
                 #print "cell ", cell.show() 
-                points += cell.retrieve()
+                points += cell.retrieve(get_vertices=get_vertices)
         return points, self.branch
 
@@ -200 +201 @@
 
 
-    def retrieve(self):
+    def retrieve_triangles(self):
+        """return a list of lists. For the inner lists,
+        The first element is the triangle index,
+        the second element is a list.for this list
+           the first element is a list of three (x, y) vertices,
+           the following elements are the three triangle normals.
+
+        This info is used in searching for a triangle that a point is in.
+
+        Post condition
+        No more points can be added to the quad tree, since the
+        points data structure is removed.
+        """
+        # FIXME Tidy up the structure that is returned.
+        # if the triangles att has been made
+        # return it.
+        if not hasattr(self,'triangles'):
+            # use a dictionary to remove duplicates
+            triangles = {}
+            verts = self.retrieve_vertices()
+            # print "verts", verts
+            for vert in verts:
+                triangle_list = self.__class__.mesh.get_triangles_and_vertices_per_node(vert)
+                for k, _ in triangle_list:
+                    if not triangles.has_key(k):
+                        # print 'k',k
+                        tri = self.__class__.mesh.get_vertex_coordinates(k)
+                        n0 = self.__class__.mesh.get_normal(k, 0)
+                        n1 = self.__class__.mesh.get_normal(k, 1)
+                        n2 = self.__class__.mesh.get_normal(k, 2) 
+                        triangles[k]=(tri, (n0, n1, n2))
+            self.triangles = triangles.items()
+            # Delete the old cell data structure to save memory
+            del self.points 
+        return self.triangles
+            
+    def retrieve_vertices(self):
+         return self.points
+
+
+    def retrieve(self, get_vertices=True):
          objects = []
          if self.children is None:
-             objects = self.points
+             if get_vertices is True:
+                 objects = self.retrieve_vertices()
+             else:
+                 objects =  self.retrieve_triangles()
          else:  
              for child in self:
-                 objects += child.retrieve()
-         return objects  
-
+                 objects += child.retrieve(get_vertices=get_vertices)
+         return objects
+        
 
     def count(self, keywords=None):

anuga_core/source/anuga/utilities/test_quad.py

@@ -1 +1 @@
 import unittest
+from Numeric import array, allclose
+
 from quad import Cell, build_quadtree
-
-#from domain import *
 from anuga.abstract_2d_finite_volumes.general_mesh import General_mesh as Mesh
 
@@ -24 +24 @@
 
         points = [a, b, c, d, e, f, g, h]
+        
         #bac, bce, ecf, dbe, daf, dae
         vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]]
@@ -56 +57 @@
         self.cell.split(4)
 
-        result =  self.cell.search(x = 1, y = 101)
+        result =  self.cell.search(x = 1, y = 101, get_vertices=True)
         assert type(result) in [types.ListType,types.TupleType],\
                                 'should be a list'
@@ -126 +127 @@
 
 
-        result = Q.search(3, 105)
+        result = Q.search(3, 105, get_vertices=True)
         assert type(result) in [types.ListType,types.TupleType],\
                                 'should be a list'
@@ -151 +152 @@
         Q = build_quadtree(mesh)
 
+    def test_retrieve_triangles(self):
+
+        cell = Cell(0, 6, 0, 6, 'cell', max_points_per_cell=4)
+
+        p0 = [2,1]
+        p1 = [4,1]
+        p2 = [4.,4]
+        p3 = [2,4]
+        p4 = [5,4]
+
+        points = [p0,p1,p2, p3, p4]
+        #
+        vertices = [[0,1,2],[0,2,3],[1,4,2]]
+
+        mesh = Mesh(points, vertices)
+
+        Q = build_quadtree(mesh)
+        results = Q.search(5,1)
+        assert len(results),2
+        #print "results", results
+        #print "results[0][0]", results[0][0]
+        assert results[0],0
+        assert results[1],2
+        assert results[0][1],[[ 2.,  1.],
+                     [ 4.,  1.],
+                     [ 4.,  4.]]
+        assert results[1][1],[[ 4.,  1.],
+                     [ 5.,  4.],
+                     [ 4.,  4.]]
+        # this is the normals
+        assert results[0][1][1],[[1.,  0.],
+                     [-0.83205029,  0.5547002],
+                     [ 0.,  -1.]]
+                     
+        # assert allclose(array(results),[[[ 2.,  1.],
+        #[ 4.,  1.], [ 4.,  4.]], [[ 4.,  1.],[ 5.,  4.],[ 4.,  4.]]] )
+        results = Q.search(5,4.)
+        ### print "results",results 
+        # results_dic={}
+        # results_dic.update(results)
+        assert len(results),3
+        #print "results_dic[0]", results_dic[0]
+        assert results[0][1],[[ 2.,  1.],
+                     [ 4.,  1.],
+                     [ 4.,  4.]]
+        assert results[1][1],[[ 2.,  1.],
+                     [ 4.,  4.],
+                     [ 2.,  4.]]
+        assert results[2][1],[[ 4.,  1.],
+                     [ 5.,  4.],
+                     [ 4.,  4.]]
+        #assert allclose(array(results),[[[ 2.,  1.],[ 4.,  1.], [ 4.,  4.]]
+         #                               ,[[ 2.,  1.],[ 4.,  4.], [ 2.,  4.]],
+        #[[ 4.,  1.],  [ 5.,  4.], [ 4.,  4.]],
+         #                               [[ 4.,  1.], [ 5.,  4.], [ 4.,  4.]]])
+        
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":
 
     mysuite = unittest.makeSuite(Test_Quad,'test')
+    # mysuite = unittest.makeSuite(Test_Quad,'test_retrieve_triangles')
     runner = unittest.TextTestRunner()
     runner.run(mysuite)