Changeset 4593

Timestamp:

Jul 5, 2007, 12:43:43 PM (17 years ago)

Author:

ole

Message:

Work on search_functions testing

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/general_mesh.py (modified) (4 diffs)
• abstract_2d_finite_volumes/test_general_mesh.py (modified) (2 diffs)
• fit_interpolate/search_functions.py (modified) (6 diffs)
• fit_interpolate/test_search_functions.py (modified) (3 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -262 +262 @@
         
 
-    def get_vertex_coordinates(self, absolute=False):
+    def get_vertex_coordinates(self,
+                               triangle_id=None,
+                               absolute=False):
         """Return vertex coordinates for all triangles. 
         
@@ -269 +271 @@
         M the number of triangles in the mesh.
 
+        if triangle_id is specified (an integer) the 3 vertex coordinates
+        for triangle_id are returned.
+        
         Boolean keyword argument absolute determines whether coordinates
         are to be made absolute by taking georeference into account
@@ -278 +283 @@
             if not self.geo_reference.is_absolute():
                 V = self.geo_reference.get_absolute(V)
+
+        if triangle_id is None:        
+            return V
+        else:
+            i = triangle_id
+            msg = 'triangle_id must be an integer'
+            assert int(i) == i, msg
+            assert 0 <= i < self.number_of_triangles
             
-        return V
+            i3 = 3*i
+            return array([V[i3,:], V[i3+1,:], V[i3+2,:]])
 
 
@@ -288 +302 @@
         """
 
-        V = self.get_vertex_coordinates(absolute=absolute)
-        return V[3*i+j, :]
+        msg = 'vertex id j must be an integer in [0,1,2]'
+        assert j in [0,1,2], msg
+        
+        V = self.get_vertex_coordinates(triangle_id=i,
+                                        absolute=absolute)
+        return V[j,:]
+    
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -41 +41 @@
                 k = triangles[i,j]  #Index of vertex j in triangle i
                 assert allclose(V[3*i+j,:], nodes[k])
+
+
+    def test_get_vertex_coordinates_triangle_id(self):
+        """test_get_vertex_coordinates_triangle_id
+        Test that vertices for one triangle can be returned.
+        """
+        from mesh_factory import rectangular
+        from Numeric import zeros, Float
+
+        #Create basic mesh
+        nodes, triangles, _ = rectangular(1, 3)
+        domain = General_mesh(nodes, triangles)
+
+
+        assert allclose(domain.get_nodes(), nodes)
+
+
+        M = domain.number_of_triangles        
+
+        for i in range(M):
+            V = domain.get_vertex_coordinates(triangle_id=i)
+            assert V.shape[0] == 3
+
+            for j in range(3):
+                k = triangles[i,j]  #Index of vertex j in triangle i
+                assert allclose(V[j,:], nodes[k])
 
 
@@ -197 +223 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_General_Mesh,'test_vertex_value_indices')
+    suite = unittest.makeSuite(Test_General_Mesh,'test')    
+    #suite = unittest.makeSuite(Test_General_Mesh,'test_vertex_value_indices')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/fit_interpolate/search_functions.py

@@ -13 +13 @@
     """
     Find the triangle (element) that the point x is in.
+
+    Inputs:
+        root: A quad tree of the vertices
+        mesh: The underlying mesh
+        x:    The point being placed
     
-    root: A quad tree of the vertices
-    Return the associated sigma and k values
-    (and if the element was found) .
+    Return:
+        element_found, sigma0, sigma1, sigma2, k
+
+        where
+        element_found: True if a triangle containing x was found
+        sigma0, sigma1, sigma2: The interpolated values
+        k: Index of triangle (if found)
+
     """
     #Find triangle containing x:
@@ -55 +65 @@
     This is called by search_tree_of_vertices once the appropriate node
     has been found from the quad tree.
+    
 
     This function is responsible for most of the compute time in
@@ -67 +78 @@
     sigma0 = -10.0
     sigma1 = -10.0
-    k = -10.0
-    #print "*$* candidate_vertices", candidate_vertices
+    k = -10
+    
     #For all vertices in same cell as point x
     for v in candidate_vertices:
@@ -79 +90 @@
         
         #for each triangle id (k) which has v as a vertex
-
         vertexlist = mesh.get_triangles_and_vertices_per_node(node=v)
         for k, _ in vertexlist:
-            #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 vertex_points of candidate triangle k
+            xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)            
 
-            #Get the three normals
-            n0 = mesh.get_normal(k, 0)
-            n1 = mesh.get_normal(k, 1)
-            n2 = mesh.get_normal(k, 2)
+            # Get the three normals (faster than using API)
+            n0 = mesh.normals[k,0:2]
+            n1 = mesh.normals[k,2:4]
+            n2 = mesh.normals[k,4:6]            
+
             
             #Compute interpolation
@@ -97 +106 @@
             sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
 
-            # Integrity check - machine precision is too hard.
-            #epsilon = get_machine_precision()            
-            # Use single precision as we used to here
+            # Integrity check - machine precision is too hard
+            # so we use hardwired single precision 
             epsilon = 1.0e-6
             msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
@@ -105 +113 @@
             assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon, msg
 
+
             #Check that this triangle contains the data point
             
             # Sigmas are allowed to get negative within
-            # machine precision on some machines (e.g nautilus)
+            # machine precision on some machines (e.g. nautilus)
             epsilon = get_machine_precision() * 2
-            if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2>= -epsilon:
+            if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
                 element_found = True
                 break
             
         if element_found is True:
-            #Don't look for any other triangle
+            # Don't look for any other triangle
             break
+        
     return element_found, sigma0, sigma1, sigma2, k
 

anuga_core/source/anuga/fit_interpolate/test_search_functions.py

@@ -3 +3 @@
 
 import unittest
-from search_functions import *
+from search_functions import search_tree_of_vertices
+from search_functions import _search_triangles_of_vertices
+
 
 from Numeric import zeros, array, allclose
@@ -9 +11 @@
 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
 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
 
@@ -49 +51 @@
         mesh.check_integrity()
 
-        #print mesh.nodes
-        #print mesh.triangles
         root = build_quadtree(mesh, max_points_per_cell = 1)
-        #print 'root', root.show()
 
-        x = [0.7, 0.7]
+        x = [0.2, 0.7]
         found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
-        #print k
-        # What is k??
+        assert k == 1 # Triangle one
         assert found is True
+
+    def test_bigger(self):
+        """test_larger mesh
+        """
+
+        points, vertices, boundary = rectangular(4, 4, 1, 1)
+        mesh = Mesh(points, vertices, boundary)
+
+        #Test that points are arranged in a counter clock wise order
+        mesh.check_integrity()
+
+        root = build_quadtree(mesh, max_points_per_cell = 4)
+
+        for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7],
+                  [0.1,0.9], [0.4,0.6], [0.9,0.1],
+                  [10, 3]]:
+                
+            found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
+
+            if k >= 0:
+                V = mesh.get_vertex_coordinates(k) # nodes for triangle k
+                assert is_inside_polygon(x, V)
+                assert found is True
+                #print k, x
+            else:
+                assert found is False                
+
+        
+
+    def test_large(self):
+        """test_larger mesh and different quad trees
+        """
+
+        points, vertices, boundary = rectangular(10, 12, 1, 1)
+        mesh = Mesh(points, vertices, boundary)
+
+        #Test that points are arranged in a counter clock wise order
+        mesh.check_integrity()
+
+        
+        for m in range(8):
+            root = build_quadtree(mesh, max_points_per_cell = m)
+            #print m, root.show()
+
+            for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7],
+                      [0.1,0.9], [0.4,0.6], [0.9,0.1],
+                      [10, 3]]:
+                
+                found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
+
+                if k >= 0:
+                    V = mesh.get_vertex_coordinates(k) # nodes for triangle k
+                    assert is_inside_polygon(x, V)
+                    assert found is True
+                else:
+                    assert found is False                
+
+                
+
+    def test_underlying_function(self):
+        """test_larger mesh and different quad trees
+        """
+
+        points, vertices, boundary = rectangular(2, 2, 1, 1)
+        mesh = Mesh(points, vertices, boundary)
+
+        root = build_quadtree(mesh, max_points_per_cell = 4)
+
+        # One point
+        x = [0.5, 0.5]
+        candidate_vertices = root.search(x[0], x[1])
+
+        #print x, candidate_vertices
+        found, sigma0, sigma1, sigma2, k = \
+               _search_triangles_of_vertices(mesh,
+                                             candidate_vertices,
+                                             x)
+
+        if k >= 0:
+            V = mesh.get_vertex_coordinates(k) # nodes for triangle k
+            assert is_inside_polygon(x, V)
+            assert found is True
+        else:
+            assert found is False                
+
+        
+
+        # More points    
+        for x in [[0.6, 0.3], [0.1, 0.2], [0.7,0.7],
+                  [0.1,0.9], [0.4,0.6], [0.9,0.1],
+                  [10, 3]]:
+                
+            candidate_vertices = root.search(x[0], x[1])
+
+            #print x, candidate_vertices
+            found, sigma0, sigma1, sigma2, k = \
+                   _search_triangles_of_vertices(mesh,
+                                                 candidate_vertices,
+                                                 x)
+            if k >= 0:
+                V = mesh.get_vertex_coordinates(k) # nodes for triangle k
+                assert is_inside_polygon(x, V)
+                assert found is True
+            else:
+                assert found is False                
+            
+
+        
+