Changeset 4599

Timestamp:

Jul 6, 2007, 12:28:53 PM (17 years ago)

Author:

ole

Message:

More work on search functions

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/general_mesh.py (modified) (3 diffs)
• abstract_2d_finite_volumes/test_general_mesh.py (modified) (3 diffs)
• fit_interpolate/search_functions.py (modified) (3 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -403 +403 @@
         Return list of triangle_ids, vertex_ids for specified node.
         If node in None or absent, this information will be returned
-        for all (full) nodes.
+        for all (full) nodes in a list L where L[v] is the triangle
+        list for node v.
         """
 
@@ -421 +422 @@
 
                 triangle_list.append( (volume_id, vertex_id) )
+
+            triangle_list = array(triangle_list)    
         else:
             # Get info for all nodes recursively.
@@ -426 +429 @@
             # working directly with the inverted triangle structure
             for i in range(self.number_of_full_nodes):
+                
                 L = self.get_triangles_and_vertices_per_node(node=i)
+
                 triangle_list.append(L)
 

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -166 +166 @@
         # One node
         L = domain.get_triangles_and_vertices_per_node(node=2)
-
         assert allclose(L[0], [0, 2])
         assert allclose(L[1], [1, 1])
@@ -173 +172 @@
         # All nodes
         ALL = domain.get_triangles_and_vertices_per_node()
-
         assert len(ALL) == 6
         for i, Lref in enumerate(ALL):
@@ -224 +222 @@
 if __name__ == "__main__":
     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

@@ -60 +60 @@
     return element_found, sigma0, sigma1, sigma2, k
 
+
 def _search_triangles_of_vertices(mesh, candidate_vertices, x):
     """Search for triangle containing x amongs candidate_vertices in mesh
@@ -82 +83 @@
     #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.
@@ -89 +91 @@
             continue
         
-        #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 k
-            xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)            
+        # 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)
 
-            # 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]            
+        # Find triangle that contains x (if any) and interpolate
+        element_found, sigma0, sigma1, sigma2, k =\
+                       find_triangle_compute_interpolation(mesh,
+                                                           triangle_list,
+                                                           x)
 
-            
-            #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)
+        if element_found is True:
+            # Don't look for any other triangle           
+            break
 
-            # 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'\
-                  %(abs(sigma0+sigma1+sigma2-1), epsilon)
-            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)
-            epsilon = get_machine_precision() * 2
-            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
-            break
         
     return element_found, sigma0, sigma1, sigma2, k
 
 
+            
+def find_triangle_compute_interpolation(mesh, triangle_list, x):
+    """Compute linear interpolation of point x and triangle k in mesh.
+    It is assumed that x belongs to triangle k.
+    """
+
+    element_found = False
+    for k, _ in triangle_list:
+        # 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)            
+        
+        
+        # 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)
+
+        # Integrity check - machine precision is too hard
+        # so we use hardwired single precision 
+        epsilon = 1.0e-6
+        delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
+        msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
+              %(delta, epsilon)
+        assert delta < 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)
+        epsilon = get_machine_precision() * 2
+        if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
+            element_found = True
+            break
+            
+    return element_found, sigma0, sigma1, sigma2, k